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Capitulo 1 
Introduccion 



Una transicion de fase se define como un cambio brusco en la estructura interna 
y las propiedades de un sistema debido a variaciones en su entorno. Este entorno 
se caracteriza generalmente por cantidades tales como temperatura, presion, cam- 
pos electromagneticos, etc. Los ejemplos mas comunes de transiciones de fase son 
la transicion de liquido a gas, de conductor normal a superconductor, o de mate- 
rial paramagnetico a ferromagnetico. El estudio de las transiciones de fase es de 
indiscutible interes tanto teorico como tecnologico. 

Los estudios teoricos microscopicos de las transiciones de fase implican el estu- 
dio de un fenomeno producido por la interaccion simultanea de un mimero enorme 
(~ 10^^) de componentes individuales. Esto forzo el desarrollo de teorias aproxima- 
das que proporcionaban soluciones exactas solo en algunos casos simplificados. Un 
ejemplo es la teoria de campo medio para transiciones de segundo orden, vease [1] 
6 [2], introducida por L. D. Landau al final de la decada de 1950. La explicacion 
mas satisfactoria de los fenomenos crfticos fue proporcionada por el Grupo de Re- 
normalizacion (GR), en primer lugar intuido por L. P. Kadanoff [3] y finalmente 
desarroUado alrededor de 1970 en los importantes articulos de K. G. Wilson [HIS], 
ver [B] para una interesante revision historica de los logros del GR. 

La transicion de fase en un sistema puede ser descrita como una discontinuidad 
en las derivadas de su energfa libre respecto a alguna de las variables termodinamicas 
y pueden ser clasificadas de acuerdo a esto, utilizando la llamada clasificacion de 
Ehrenfest. Si la discontinuidad se presenta en la primera derivada, se denomina 
transicion de fase de primer orden, mientras que si es en la segunda derivada, se 
denomina transicion de fase de segundo orden. 

De forma general, las transiciones de fase de primer orden son casi siempre las 
que involucran un calor latente. Durante dichas transiciones, el sistema absorbe o 
libera una cantidad fija (y por lo general grande) de energia. Durante el proceso, 
la temperatura del sistema permanece constante conforme se absorbe o se libera 
calor. Ademas, las transiciones de primer orden estan asociadas a regimenes mixtos 
en los que algunas partes del sistema ban completado la transicion, mientras que 
otras no. Un ejemplo tfpico de este fenomeno es la coexistencia del regimen de baja 
temperatura del agua (hielo) y el de alta temperatura (agua liquida); el agua y el 
hielo pueden coexistir (existen los icebergs). 
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Las transiciones de fase de segundo orden son continuas en la primera derivada 
de la energi'a libre, pero presentan discontinuidades en su segunda derivada. Estas 
incluyen la transicion a la fase ferromagnetica en materiales como el liierro, donde la 
magnet izacion, que es la primera derivada de la energi'a libre con respecto a la fuerza 
del campo magnetico aplicado, aumenta de forma continua desde cero conforme la 
temperatura desciende por debajo de la temperatura de Curie. La susceptibilidad 
magnetica, la segunda derivada de la energi'a libre respecto al campo, diverge. Este 
tipo de transiciones no tiene calor latente asociado pero presenta longitudes de 
correlacion infinita. Ejemplos ti'picos de transiciones de fase de segundo orden son las 
transiciones paramagnetica-ferromagnetica y conductor-superconductor. Este tipo 
de transiciones tambien se caracteriza por comportamientos en forma de leyes de 
potencia en el punto de transicion (tambien llamado punto critico) con exponentes 
no enteros, llamados exponentes critico^ Los exponentes cn'ticos estan relacionados 
entre si por relaciones de hiperescalado - conociendo dos de los exponentes, los otros 
pueden ser deducidos. Sistemas muy diferentes pueden compartir exactamente el 
mismo conjunto de exponentes cn'ticos y se dice entonces que pertenecen a la misma 
Clase de Universalidad (CU). La CU de un sistema se define por propiedades muy 
generales tales como la simetria de la interaccion microscopica, la dimensionalidad 
espacial, o la dimensionalidad del parametro de orden, ver [8] para una revision 
exhaustiva de las CU mas habituales. 

La existencia de un parametro de orden es tambien caracteristica de las tran- 
siciones de fase. Este puede ser definido como una cantidad que es nula en una de 
las fases y no nula en la otra. Refieja el proceso de ruptura de simetria que normal- 
mente tiene lugar a traves del punto de transicion. Por ejemplo, para la transicion 
paramagneto-ferromagneto un parametro de orden valido es la magnetizacion neta 
(cero en la fase de alta temperatura y no cero en la de baja temperatura), mientras 
que para la transicion li'quido-gas es la diferencia de densidad de los dos regi'menes 
que coexisten. Otros tipos de transiciones de fase deben ser descritos por parametros 
de orden mas complejos. 

Simulaciones de Monte Carlo (MC) ban resultado ser muy utiles en esta rama 
de la Mecanica Estadfstica, ver [9] para una revision de los metodos mas populares. 
Con ellos, se puede simular la evolucion temporal de cada constituyente del sistema 
para un determinado Hamiltoniano. En nuestro caso, los sistemas se definen en 
redes de dimensionalidad espacial D, con tamaiio lineal L y condiciones de contorno 
periodicas. En cada nodo de la red se define una variable, Uamada espin, que toma 
ciertos valores (dependientes del modelo) que evolucionan con el tiempo. 

Un metodo de MC puede actualizar o bien un solo espm por iteracion, como es 
el caso del algoritmo de Metropolis o del de baiio termico [9], o bien un grupo de 
espines, como en el caso de los algoritmos de Wolff [TO] o de Swendsen-Wang [TT] . 
Estos ultimos son llamados metodos de clusters. En las proximidades del punto cn'ti- 
co, se produce el llamado Decaimiento Critico, ver por ejemplo [9j. Los tiempos de 

^ Seguiremos la nomenclatura habitual (ver por ejemplo [7]) para los exponentes cn'ticos: v 
es el exponente de la longitud de correlacion, a es el del calor especi'fico, /? el del parametro de 
orden, mientras que w es el exponente (universal) de correcciones de escala de orden dominante. 
Un exponente ligeramente distinto, la dimension anomala 77, se define en Ec. (|2.72p . 
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relajacion del parametro de orden divergen como una potencia de la longitud de co- 
rrelacion, r ~ siendo z el denominado exponente critico dindmico. Esto implica 
que el tiempo necesario para producir configuraciones estadisticamente independien- 
tes diverge cerca del punto critico para un sistema finito como r ~ L^. Los metodos 
de MC de actualizacion de un linico espfn tienen un exponente z '>2. Por lo tanto 
es muy complicado obtener datos de alta precision muy cerca del punto critico en 
sistemas grandes. Por el contrario, se obtiene un comportamiento dinamico mucho 
mejor utilizando metodos de MC de actualizacion de clusters. Con estos liltimos, 
dependiendo del modelo y de la dimensionalidad, se obtienen valores de z entre y 
1 [12]. En este trabajo hemos utilizado casi siempre algoritmos de clusters. 

Los metodos de actualizacion de espines incluyen simulaciones dentro del co- 
lectivo canonico (a una temperatura fija) y dentro del colectivo microcanonico (a 
energia fija). Respecto a las simulaciones dentro del colectivo microcanonico, hemos 
utilizado un metodo de simulacion microcanonico propuesto recientemente que per- 
mite la simulacion de sistemas de un tamario nunca antes alcanzado que realizan 
transiciones de fase de primer orden [T3] . 

Incluso con los recursos de computacion de hoy en di'a, estamos restringidos a 
simular sistemas con mas de 10^^ ordenes de magnitud menos componentes que los 
presentes en un sistema macroscopico real (con ~ 10^^ particulas). Lo unico que 
podemos hacer es simular sistemas con diferentes tamafios y tratar de extrapolar 
nuestros resultados al Limite Termodinamico {L — )■ oo). El estudio del comporta- 
miento de escala de los diferentes observables con el tamano del sistema es Uamado 
Finite-Size Scaling (FSS) y es fundamental para el estudio de las transiciones de 
fase, vease por ejemplo |7]. En este trabajo hemos usado continuamente tecnicas 
de FFS, ademas hemos realizado un estudio novedoso del FSS dentro del colectivo 
microcanonico, vease el Capitulo[2J 

Nuestro objetivo principal es el estudio de los efectos del desorden sobre las tran- 
siciones de fase. En concreto, deseamos estudiar el efecto de las impurezas congeladas 
en la transicion de material paramagnetico a ferromagnetico. La presencia de des- 
orden aleatorio en un sistema produce muchos fenomenos interesantes y fisicamente 
relevantes, lo que ha motivado extensos estudios teoricos y experimentales. Los tipos 
mas caracteristicos de sistema con desorden aleatorio son: vidrios de espin p^lTS]. 
sistemas aleatoriamente anisotropicos [TBHTS] . sistemas diluidos P^VET] y sistemas 
con los campos aleatorios [H]. En todos estos casos, existen variables aleatorias que 
caracterizan el comportamiento del sistema. 

Al modelar un sistema con desorden aleatorio se pueden utilizar dos enfoques 
diferentes. Por un lado, se puede considerar que las variables aleatorias estan en 
equilibrio termodinamico las otras variables dinamicas del sistema. Por lo que las 
variables aleatorias tambien seran "dinamicas" . Este es el Uamado desorden annea- 
led y debe ser la eleccion si deseamos modelar un sistema en el que los tiempos 
caracteristicos de la dinamica del desorden sean comparables con los tiempos carac- 
teristicos de la dinamica de las variables originales, como seria el caso por ejemplo de 
una disolucion de dos liquidos. Por otro lado, se puede considerar que las variables 
aleatorias no evolucionan en el tiempo, sino que estan congeladas. Es el Uamado des- 
orden congelado. Esta ultima es una alternativa perfectamente valida por ejemplo si 
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queremos modelar sistemas magneticos con impurezas. En este caso, el comporta- 
miento magnetico se debe a los espines de los electrones no apareados en las capas 
atomicas exteriores, mientras que las impurezas son atomos sin electrones no apa- 
reados. Se sabe que la dinamica de los electrones es ordenes de magnitud mas veloz 
que la dinamica de los nucleos de modo que se puede considerar a los atomos de im- 
purezas congelados en el tiempo. En Ref. [22] se presenta un analisis mas detallado 
de esta cuestion. Cuando se considera el desorden congelado se generan diferentes 
configuraciones espaciales aleatorias del desorden (llamadas muestras). Los espines 
de cada muestra evolucionaran independientemente mientras que las impurezas per- 
manecen fijas. Para extraer la informacion de un determinado observable, en primer 
lugar realizamos el promedio de su evolucion temporal en el interior de cada muestra 
(en lo sucesivo denominado promedio termal y denotado por brackets) y luego rea- 
lizamos el promedio de lo anterior entre todas las muestras (denominado promedio 
muestral y denotado por un suprarrayado), el promedio doble es entonces denotado 
por (■■■). 

Uno de los resultados de mayor importancia en el estudio de sistemas desordena- 
dos es el criterio de Harris [22] , vease el Apendice El El criterio sefiala que si el calor 
especi'fico en el sistema puro diverge (el exponente critico ctpum es mayor que cero), 
el desorden cambiara el comportamiento critico del modelo, es decir, aparecera una 
nueva CU. En este caso, se dice que el desorden es relevante. Por el contrario, si el 
calor especi'fico no diverge en el sistema puro (ttpuro < 0) los exponentes criticos del 
sistema desordenado no cambiaran. Es este caso, se dice que el desorden es irrele- 
vante. En el presente trabajo recomprobaremos la validez del criterio para el modelo 
de Heisenberg tridimensional con dilucion por sitios. 

Otra cuestion muy interesante que surge cuando se estudian sistemas diluidos 
es la cuestion del autopromedio. El valor medio de un observable O en una red de 
tamano lineal L es diferente para cada realizacion del desorden (en nuestro caso, 
para cada distribucion espacial de los sitios no magneticos), por lo que es una varia- 
ble estocastica caracterizada por un promedio sobre el desorden O y una varianza 
(Z\0)2 = O^-O . Se dice que un sistema exhibe autopromedio para el observable O 
si AO/O tiende a cero cuando L — )■ oo. Cuando un sistema diluido no autopromedia 
los estudios numericos se hacen muy dificiles: incluso fijando la temperatura crftica 
al valor correcto para L — )■ oo, hacer el sistema mas grande no proporciona una gran 
mejora estadistica. El autopromedio de las propiedades de los sistemas desordena- 
dos genera gran interes, reflejado en numerosos trabajos tanto numericos [2H426] 
como analiticos [27H29] . En esta tesis se estudiara el autopromedio tanto de la sus- 
ceptibilidad del modelo de Heisenberg tridimensional, como del calor latente y la 
tension superficial del modelo de Potts tridimensional, ambos modelos con dilucion 
por sitios. 

Hemos estudiado numericamente varios modelos con desorden aleatorio que pre- 
sentaban importantes cuestiones abiertas. Nuestra colaboracion ha producido las 
publicaciones recogidas en [2DH32] aunque en el presente trabajo solo se presentan 
los resultados de las referencias [SUHSS] . 

La disposicion del resto de esta tesis doctoral es la siguiente. En el Capitulo [2] 
estudiamos las propiedades de escala, tanto del modelo de Potts puro con cuatro 
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estados [Q = 4) en D = 2 como del modelo de Ising puro en D = 3. Estos modelos 
realizan transiciones de fase de segundo orden bien conocidas con calores especifi- 
cos divergentes. Hemos simulado ambos modelos utilizando el metodo de simulacion 
microcanonico presentado en [13]. Obtuvimos fuertes evidencias de la bondad de 
dicho metodo a traves de la comparacion con los resultados mas recientes [8lll3]. El 
Capitulo 13] lo dedicamos al estudio de los efectos de la dilucion en un sistema que 
realiza una transicion de fase de primer orden fuerte: el modelo de Potts tridimen- 
sional con Q = 4 J Q = 8 estados. Utilizando el metodo de simulacion comprobado 
en el Capitulo [2], fuimos capaces de simular sistemas con mas de 10^ componentes, 
multiplicando por un factor de 100 el numero de componentes de los trabajos mas 
recientes [HH15] . En el Capitulo |l]estudiamos las propiedades de autopromedio del 
modelo de Heisenberg tridimensional diluido por sitios, donde existen dos escenarios 
en conflicto que afirman que la susceptibilidad es [27] o no es pHj autopromedian- 
te. Ademas obtendremos informacion acerca de la validez del criterio de Harris. En 
el Capitulo [5] exponemos nueva informacion acerca de los exponentes criticos de 
los terminos logaritmicos del modelo de Ising en cuatro dimensiones con dilucion 
por sitios. Hemos logrado discriminar entre cinco diferentes teorias [17H5T] median- 
te simulaciones numericas de gran envergadura. Presentamos nuestras conclusiones 
generales en el Capitulo (6] Tambien se exponen varios apendices tratando de am- 
pliar la informacion sobre algunas de las herramientas utilizadas mas importantes 
o innovadoras. En el Apendice \M explicamos con detalle el criterio de Harris. En 
el Apendice [B] presentamos brevemente las tecnicas populares del FSS, asi como 
el Metodo de los Cocientes, una tecnica que permite el calculo de los exponente 
criticos partiendo de los datos obtenidos en sistemas finitos. En el Apendice O se 
describen dos cuestiones muy importantes para la simulacion de sistemas dinamicos 
con metodos de MC: los tiempos de autocorrelacion y la estimacion de errores. El 
Apendice [D] esta dedicado a describir los diferentes metodos para la extrapolacion 
en temperatura de los resultados obtenidos en una simulacion canonica mientras 
que en el Apendice [E] se describe la obtencion de la construccion de Maxwell, muy 
litil para el estudio de transiciones de fase de primer orden. En el Apendice |F] se 
describe el enfoque introducido por Lee y Yang para describir las transiciones de 
fase, formulando su importante teorema. Tambien discutimos en este apendice de 
la distribucion de los ceros de Lee- Yang sobre el circulo unidad. Por ultimo, dedi- 
camos el Apendice [G] a la descripcion de la infraestructura de supercomputacion 
IBERCIVIS, que ha sido crucial para completar algunas partes de este trabajo. 
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Chapter 1 
Introduction 



A phase transition is defined as a sharp change in the internal structure and prop- 
erties of a system due to variations in its environment. This environment is usually 
characterised by quantities such as temperature, pressure, electromagnetic fields, 
etc. Common examples of phase transitions are the transition from liquid to gas, 
from normal conductor to superconductor, and from paramagnet to ferromagnet. 
The study of phase transitions is of major interest both theoretically and technolo- 
gically. 

Theoretical microscopic studies of phase transitions involve the study of a phe- 
nomenon produced by the simultaneous interaction of an enormous number (~ 10^^) 
of individual components. This forced the development of approximate theories that 
produce exact solutions only for some simplified cases. An example is the mean-field 
theory for second-order phase transitions, see for instance PQ or introduced by 
L. D. Landau in the late 1950s. The most satisfying explanation of critical phe- 
nomena was provided by the Renormalization Group (RG) picture, first intuited 
by L. P. Kadanoff [3] and finally developed around 1970 in the landmark papers of 
K. G. Wilson ^E], see [6] for an interesting historical review of RG achievements. 

A phase transition in a system can be described as a discontinuity in the de- 
rivatives of its free energy with respect to some thermodynamic variable, and can 
be classified according to this by using the so-called Ehrenfest classification. If the 
discontinuity is in the first derivative it is called a first-order phase transition, while 
if it is in the second derivative it is called a second-order phase transition. 

More generally, first-order phase transitions are usually those involving a lat- 
ent heat. During such a transition, the system either absorbs or releases a fixed 
(and typically large) amount of energy. During the process, the temperature of the 
system remains constant as heat is added or released. In addition, first-order trans- 
itions are associated with mixed-phase regimes in which some parts of the system 
have completed the transition while others have not. A typical example of this phe- 
nomenon is the coexistence of the low temperature regime of water (ice) and the 
high temperature one (liquid water); water and ice can and do coexist (there exist 
icebergs) . 

Second-order phase transitions are continuous in the first derivative but exhibit 
discontinuities in a second derivative of the free energy. These include the ferro- 
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magnetic phase transition in materials such as iron, where the magnetisation, which 
is the first derivative of the free energy with respect to the apphed magnetic field 
strength, increases continuously from zero as the temperature is lowered below the 
Curie temperature. The magnetic susceptibility, the second derivative of the free 
energy with respect to the field, diverges. These have no associated latent heats, 
but present infinite correlation lengths. Examples of second-order phase transitions 
are the paramagnetic-ferromagnetic and the conductor-superconductor transitions. 
This kind of transition is also characterised by power-law behaviour at the trans- 
ition point (also-called critical point) with non-integer exponents, thus called critical 
exponent^ The critical exponents are related to each other by hyperscaling rela- 
tions - knowing two of the exponents the others can be deduced. Quite different 
systems can share exactly the same set of critical exponents and are said to belong 
to the same Universality Class (UC). The UC of a system is defined by very general 
properties such as symmetry of the microscopic interaction, dimensionality of the 
space, or dimensionality of the order parameter, see [8] for an exhaustive review of 
the most usual UC's. 

The existence of an order parameter is also characteristic of phase transitions. 
It can be defined as a quantity that is null in one of the phases and non-null in the 
other, and reflects the symmetry-breaking process that usually takes place across 
the transition point. For example, for the ferromagnetic-paramagnetic transition a 
valid order parameter is the net magnetisation (zero in the high temperature phase 
and non-zero in the low-temperature one), while for the liquid-gas transition it is 
the density difference of the two co-existing regimes. Other kinds of phase transition 
must be described by more complex order parameters. 

Monte Carlo (MC) simulations have proved very useful in this branch of Stat- 
istical Mechanics, see |9] for a review of the most popular methods. With them, 
one can simulate the evolution in time of each constituent of the system for a given 
Hamiltonian. In our case the systems are defined on lattices of spatial dimensionality 
D, with linear size L and periodic boundary conditions. On each node of the lattice 
we define a variable, called spin, that takes on values (depending on the model) that 
evolve in time. 

An MC spin update method can change either just one spin per iteration, as 
is the case of the Metropolis or heat-bath algorithms [9], or a cluster of spins, as 
is the case of the Wolff [10] and Swendsen-Wang |T1] algorithms. The latter are 
thus called cluster methods. At the critical point is presented the so-called Critical 
Slowing Down, see for example [9]. The relaxation time of the order parameter 
diverges as a power of the correlation length, r ~ with z being the dynamic 
critical exponent. This roughly implies that the time needed to produce statistically 
independent configurations diverges at the critical point for a finite system as r ~ L^. 
Single-spin MC update methods have an exponent z >2. Therefore is very hard to 
obtain high-precision data very close to the critical point on large systems. However, 

^We follow the standard terminology (see e.g [7]) for the critical exponents: z/ is the exponent 
for the correlation length, a that of the specific heat, /3 that of the order parameter, while uj 
is the (universal) leading-order scaling-corrections exponent. A slightly different exponent, the 
anomalous dimension rj, is defined in Eq. (j2.72p . 



16 



cluster MC update methods produce a much better dynamic behaviour. Depending 
on the model and dimensionality, cluster methods have z values between and 1 [H] . 
In this work we have used basically cluster methods. 

Spin update methods include simulations within the canonical ensemble (at fixed 
temperatures) and within the microcanonical ensemble (at fixed energies). With 
respect to simulations in the microcanonical ensemble, we have exploited a recently 
proposed microcanonical simulation method that allows the simulation of systems 
with a size never reached before that undergo first-order phase transitions [13]. 

Even with today's computing resources, we are restricted to simulating systems 
with more than 10^^ orders of magnitude fewer components than the real mac- 
roscopic system (with ~ 10^^ particles). The only thing we can do is to simulate 
systems with different sizes and try to extrapolate the results to the Thermodynamic 
Limit [L — )■ oo). The study of the scaling behaviour of the different observables with 
system size is called Finite-Size Scaling (FSS) and is fundamental for the study of 
phase transitions, see for example In this work we have continually used FSS 
techniques, as well as performing a novel study of FSS within the microcanonical 
ensemble, see Chap. [H 

Our main objective is the study of the effects of disorder in phase transitions. In 
particular, we study the effect of quenched impurities on the paramagnet-ferromagnet 
transition. The presence of random disorder in a system produces many interesting 
and physically relevant phenomena which have motivated extensive theoretical and 
experimental studies. The most typical types of system with random disorder are: 
spin glasses [IHIIS], random anisotropic systems [TBHTS] . dilute systems [TUIET] . 
and systems with random fields [14J. In all these cases there exist random variables 
characterising the behaviour of the system. 

When modelling a randomly disordered system one can use either of two ap- 
proaches. On the one hand, one can consider the random variables in thermody- 
namic equilibrium with the other dynamic variables of the system. Thus the random 
variables will also be "dynamic" . This is the so-called annealed disorder and should 
be the choice if we model a system in which the characteristic times for the dynamics 
of the disorder are comparable with the characteristic time of the original dynamic 
variables, as would be the case for example of a solution of two liquids. On the 
other hand, one can consider that the random variables do not evolve in time, but 
are frozen. This is the so-called quenched disorder. The latter is a perfectly valid 
alternative for example if we want to model magnetic systems with impurities. In 
this case the magnetic behaviour is due to the spins of the unpaired electrons in the 
outer atomic shells while the impurities are whole atoms with no unpaired electrons. 
It is known that the dynamics of the electrons is orders of magnitude faster than 
the dynamics of the nuclei, so that we can perfectly consider the impurity atoms 
as frozen in time. In Ref. [22] there is presented a more detailed discussion of this 
issue. When considering quenched disorder we will generate different random spatial 
configurations of the disorder (called samples). Within each sample the spins evolve 
independently but the disorder is fixed. To extract information of a given observ- 
able, first we perform the average of its temporal evolution within each sample (in 
the following called thermal average and denoted by angle brackets) and afterwards 
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we perform the sample average (denoted by an overline), the double average is then 
denoted by (■■■). 

One of the main results in disordered systems is the Harris criterion [23], see 
Appendix |X1 It states that if the specific heat diverges in the pure system (the 
critical exponent, Opure, is greater than zero), the disorder will change the critical 
behaviour of the model, i.e a new UC will appear. In this case it is said that the 
disorder is relevant. Conversely, if the specific heat does not diverge in the pure 
system (opure < 0) the critical exponents of the disordered system will not change. 
In this case it is said that the disorder is irrelevant. We will recheck in the present 
work the validity of the criterion for the three-dimensional site-diluted Heisenberg 
model. 

Another very interesting question arising when studying dilute systems is the 
issue of self- averaging. The mean value of a quantity (9 on a lattice of linear size L is 
different for each realization of the disorder (in our case, for each spatial distribution 
of the non-magnetic sites). Therefore it is a stochastic variable characterised by an 
average over the disorder O and a variance {AOy = O'^ — O . It is said that a 
system is self-averaging for the quantity O if AO/O goes to zero when L — )■ oo. 
When a dilute system is not self- averaging, numerical studies become very difficult: 
even fixing the critical temperature to the correct value for L — > oo making the 
system larger does not much improve the statistics. The self-averaging properties 
of disordered systems have generated much interest, refiected in numerous works 
both numerical and analytical [271 - I29] . In this work we will study the self- 

averaging properties both of the susceptibility of the three-dimensional site-diluted 
Heisenberg model and of the latent heat and surface tension of the three-dimensional 
site-diluted Potts model. 

We have numerically studied randomly disordered models presenting important 
open issues. Our collaboration has produced the papers of Refs. although 
in the present work we only present the results of Refs. [SOHSSI- 

The organisation of the rest of this PhD thesis is as follows. In Chapter [2] we 
study the scaling properties both of the four-state {Q = 4) pure Potts model in 
D = 2 and of the pure Ising model in D = 3. These models undergo well-known 
second-order phase transitions with diverging specific heats. We have simulated 
them using the microcanonical simulation method presented in [13] obtaining strong 
evidence for the goodness of our approach by comparing it with the most recent 
results [8E3] . Chapter [3]is devoted to the study of the effects of dilution on a system 
performing a strong first-order phase transition: the three-dimensional Potts model 
with Q = 4 and Q = 8 states. Using the simulation method studied in Chapter [21 
we will be able to simulate systems with more than 10^ components, multiplying by 
a factor of 100 the number of components reached in the most recent work [Hl - HG] . 
In Chapter H] we study the self-averaging properties of the three-dimensional site- 
diluted Heisenberg model, where there exist two confiicting results stating that the 
susceptibility is [27] or is not [28] a self-averaging quantity. We will also obtain 
information about the validity of the Harris criterion. In Chapter [5] we report 
novel information about the critical exponents of the logarithmic terms of the four- 
dimensional site-diluted Ising model, we try to discriminate between five different 
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theories |17H^ by using high-statistics MC simulations. We present our conclusions 
in Chapter O We also present some appendices with the aim of extending some of 
the most important or innovative tools used. In Appendix |A] we explain in detail the 
Harris criterion. In Appendix |B] we briefly present the popular FSS techniques and 
the Quotient Method, a technique that allows the computation of critical exponents 
from data obtained in finite systems. In Appendix O we describe two important 
issues when simulating dynamical systems with MC methods - autocorrelation times 
and error estimates. Appendix [D] is devoted to describing the different methods to 
temperature-extrapolate the results obtained in a canonical MC simulation, and 
Appendix [E] describes the derivation of the Maxwell construction, which is very 
useful in the study of first-order phase transitions. In Appendix [F] we describe 
the approach introduced by Lee and Yang to describe phase transitions, formulating 
their landmark theorem. We also discuss in this appendix the distribution of the LY- 
zeros on the unit circle. Finally Appendix [G] describes the IBERCIVIS computing 
infrastructure, which has been crucial for the completion of some parts of this work. 
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Chapter 2 

Microcanonical Finite-Size Scaling 



2.1 Introduction 

The canonical ensemble enjoys a predominating position in Theoretical Physics due 
to its many technical advantages (convex effective potential in finite systems, eas- 
ily derived fluctuation-dissipation theorems, etc.). This somewhat arbitrary choice 
of ensemble is justified by the ensemble equivalence property, which holds in the 
Thermodynamic Limit (TL) for systems with short-range interactions. 

However, in spite of this long-standing bias in favour of the canonical ensemble, 
the canonical analysis of phase transitions is not simpler. The advantages of mi- 
crocanonical analyses of first-order phase transitions have long been known [T3l[52] , 
and indeed become overwhelming in the study of disordered systems [31]. Further- 
more, the current interest in mesoscopic or even nanoscopic systems, where ensemble 
equivalence does not hold, provides ample motivation to study other statistical en- 
sembles, in particular the microcanonical ensemble [53]. Besides, microcanonical 
Monte Carlo [51] is now as simple and efficient as its canonical counterpart (even 
microcanonical cluster algorithms are available [13] )• Under such circumstances, it 
is of interest to extend Finite-Size Scaling (FSS) [Tt lS^HFT] to the microcanonical 
framework for systems undergoing a continuous phase transition. 

The relation between the microcanonical and the canonical critical behaviour is 
well understood only in the TL. A global constraint modifies the critical exponents, 
but only if the specific heat of the unconstrained system diverges with a positive 
critical exponent a > [58] (however, see [59]). This fact is explained in detail 
in Sec. 12.2.11 The modification of the critical exponents, termed Fisher renormal- 
ization, is very simple. Let L be the system size, and consider an observable O 
(for instance, the susceptibility) whose diverging behaviour in the infinite-volume 
canonical system is governed by the critical exponent xq 

(0)i™^'«W°' t = ^^^. (2.1) 

c 

Now, let e be the internal energy density and Cc = (e)^^"J|^rf • Consider the microca- 
nonical expectation value of the same observable O in Eq. (12. ip . but now at fixed 
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energy e. The scaling behaviour f|2.ip translates to 

(0)^=00.6 oc le-ecT'^O''", xo,m = -^^- (2.2) 

1 — a 

We will denote the microcanonical exponents with the subindex "m". Hence, the 
Fisher renormalization of the correlation length exponent u, is u — )■ = z^/(l — a), 
that of the order parameter exponent is (3 — )■ = — a), etc. On the other 
hand, the anomalous dimension, defined in Eq. (12.721) . is invariant under Fisher 
renormalization |58j, i.e. rj = rj^. See also [60] for a recent extension of Fisher 
renormalization to the case of logarithmic scaling corrections. 

As for systems of finite size, the microcanonical FSS [611463] is at the level of 
an ansatz. This ansatz is obtained from the canonical one merely by replacing the 
free-energy density by the entropy density, and using Fisher renormalised critical 
exponents. The microcanonical ansatz reproduces the canonical one [64j, and has 
been the subject of some numerical testing [63(165] . Furthermore, systems undergoing 
Fisher renormalization (due to some global constraint other than the energy) do seem 
to obey FSS as well [66j. 

A difficulty lies in the fact that the current forms of the microcanonical FSS 
ansatz (FSSA) [6TH63] are in a somewhat old-fashioned form. Indeed, they are 
formulated in terms of quantities such as Cc or the critical exponents, which are 
not accessible in the absence of an analytical solution. In this respect, a great 
step forward was achieved in a canonical context [67| when it was realized that the 
finite-lattice correlation length [68] allows one to formulate the FSSA in terms of 
quantities computable in a finite-lattice. This formulation made it practical to ex- 
tend Nightingale's phenomenological renormalization ^69j to space dimensions D > 2 
(the so-called quotient method [70]). 

Here, we will extend the microcanonical FSSA to a modern form, allowing us 
to use the quotient method. We will test numerically this extended FSSA in two 
models with a > 0, hence undergoing non-trivial Fisher renormalization, namely 



^In the particular case of the fixed-energy constraint, Eq. (|2.2p follows from (j2.ip and from the 
ensemble equivalence property 

/ canonical ^ /^xcanonical 

Indeed, it suffices to notice that {C{T) is the canonical specific heat, C oc |i|~"), 



/ dTC(r) cx =^ |t| (X |e-ec|T^ 

Jt,. 



The only exponent whose renormalization is not clear at this point, is a itself, for the energy is 
not a dynamical variable but a parameter in this ensemble. If one chooses to define am as the 
critical exponent corresponding to dt/de, the correspondence with Fisher renormalization, am = 
— a/(l — a) becomes complete. In fact, see concluding paragraph in Sect. 12.2.21 the microcanonical 
dt/de behaves as the canonical l/[de/dt] . Note that in the above expressions we disregarded 
subdominant terms such as the contribution of the analytical background in the specific heat. 
Such terms are subdominant only if a > 0. In case a were negative, the asymptotic dominance is 
different. The specific heat at Tc is dominated by the analytical background. As a consequence 
\t\ ~ |e — ed and none of the exponents (not even a) gets renormalised. 
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the D = 3 ferromagnetic Ising model, and the D = 2 four-state ferromagnetic Potts 
model. The Potts model has the added interest of undergoing, in its canonical form, 
quite strong logarithmic corrections to scaling that are nevertheless under relatively 
strong analytical control |13]. It will therefore be quite a challenge to control the 
logarithmic corrections in the microcanonical setting. 

2.2 Analytical Framework 

2.2.1 Fisher Renormalization of Critical Exponents 

In 1976, an important paper of M. E. Fisher [52] established a set of elegant rela- 
tionships describing the effects of constrained hidden variables on the critical ex- 
ponents. The original theory was developed to explain the significant deviations of 
the theoretical predictions (basically from the Ising model) from the experimental 
measurements of critical exponents. These deviations were attributed to some extra 
"hidden" degrees of freedom, present in the real system but not in the oversimplified 
ideal theoretical model. Models (like Ising or Potts) are somewhat gross idealisa- 
tions of real fluids or magnets, and can be said to lack sufficient internal degrees 
of freedom. In addition, some experiments are unavoidably different from the ideal 
system, for example due to the presence of defects or impurities (such as quenched 
magnetic impurities or non-uniform isotopic composition). 

Apart from the exactly soluble models introduced in the original work [58j, the 
formalism introduced by Fisher has provided explanations of numerous phenomena 
and behaviours, both theoretical and experimental. To cite some examples, there 
are studies of the phase transition of constrained uniaxial dipolar ferromagnets [7T] 
and of the random-field antiferromagnet with competing interactions [72], efforts 
made to distinguish between the Random Field Ising Model (RFIM) and the Dilute 
Antiferromagnetic Model (DAFM) under an applied field [73] , and the study of the 
tricritical point of the Blume-Capel model in three dimensions [71] related with the 
superfiuid A transition in ^He-^He mixtures in confined films [7S]. We would also 
point to the agreement of Fisher theory with the results for compressible systems, 
theoretically for both the Ising [77] and the 0^ [TS] models and experimentally for 
ammonium chloride at high pressures [79], see also [76] for a study of the tricritical 
point in compressible systems. 

The situation was described as follows: firstly, there is an "ideal" system with 
known variables characterised by the ideal critical exponents a,/3,7, . . .; secondly, 
the "real" system has some "hidden" degrees of freedom which fiuctuate but remain 
in equilibrium with the known variables; finally, the hidden variables are subject to 
some form of constraint (for example, the total number of impurity atoms must re- 
main fixed). The critical exponents of the real system are denoted by ax, Pxylx, ■ ■ ■ 

In the following we will describe the relationships between the two sets of critical 
exponents (now called Fisher renormalization), following in some points the recent 
work of Kenna et al. [60] in which there also can be found the corresponding set of re- 
lationships for the logarithmic correction exponents. We will focus the discussion on 
describing the temperature transition in a ferromagnet at the Curie point, although 
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the results arc perfectly valid for other kinds of phase transitions (antiferromagnets, 
gas-liquid transitions, binary fluids, etc.). 

The starting point is the description of the ideal system in terms of the free 
energy: 

f-Mt,h), (2.3) 

where t = {T — Tc)/Tc is the reduced temperature and h is the "field". The field 
is related in general with the order parameter, a, describing the transition under 
study through 

= = -(I).- P-^) 

that in the ideal case becomes 

. = = -(f)^. (2.5) 

The order parameter can be the magnetisation in a ferromagnet, the sublattice 
magnetisation in an antiferromagnet, the density in a gas-liquid transition, etc. The 
ideal system will undergo a phase transition at T°. Thus below this temperature 
there will be a non- vanishing order parameter (spontaneous magnetisation) 

Aao = lim hao{t, h) - ao{t, -h)] , (2.6) 

which vanishes at the critical point as: 

Aao ~ l^l^^ , {t 0") . (2.7) 
The second field derivative is also described by its corresponding critical exponent 

. (*-0). (2.8) 

Finally the critical behaviour of the specific heat 

co(t)~|tr" , (t^O), (2.9) 
means that, in the absence of a field, the ideal free energy behaves for small t as 

fo{t, 0) = A± + A±|t| + A2±\t\' + B±\t\^-" + 0(|t|-^) , (2.10) 

where the ± depends on the sign of t. 

We assume that the real system is derived from the ideal system by the intro- 
duction of a new "hidden" thermodynamic variable x which is the conjugate of a 
force u, such that the thermodynamic potential becomes / = f{t, h, u), with 
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The constraint in the hidden variable is written as 

x{t,h,u) = X{t,h,u) , (2.12) 

where X{t, h,u) is assumed to be an analytical function. Let us introduce the basic 
hypothesis (suggested in part by the analytical results of some soluble models [58]) 
that the free energy of the constrained system can be written in terms of the ideal 
free energy /q as 

fit, h, u) = fo{t*{t, h, u), h*{t, h, u)) + g{t, h, u) , (2.13) 

where t*, h*, and g are analytic functions of their arguments. I.e., we assume that 
the total free energy consists of a "regular background" contribution g{t,h,u) plus 
a "singular contribution" derived from the ideal free energy /o(t, /i) by a smooth 
transformation of the temperature, t, and the field, h, to the modified versions t* 
and h*. Furthermore, the transition must remain ideal if observed at fixed force u, 
and the ideal free energy fo(t, h) is recovered when m = 0. 

To simplify the discussion we assume that the hidden degrees of freedom are 
neutral in the sense that they do not bias the value of the external field at the 
transition. I.e., the transition still occurs at /i = and 

h*{t,h,u) = hj{t,h,u) . (2.14) 

With the previous assumptions, one can obtain from Eqs. (12.61] , (12.131) . and (12.141) 
that the order parameter of the constrained system behaves as 

Aa = Aao{t*{t, 0, u))J{t, 0, u) , (2.15) 

where we have made use of the continuity of all the analytical functions for /i — > 
including the internal energy of the ideal system, cq = dfo/dt. 

We can calculate the internal energy for zero field using Eq. (12.101) as 

eo(t, 0) = ^^^^ = Ao + A\t\+B\t\'-'^ + ---, (2.16) 

with Aq = ±^1-1-, Ai = 2A2± and B = ±{2 — a)B± and where ■ ■ ■ represents 
higher-order terms. 

We can also obtain from Eqs. (12.111) and (12.131) for — )■ 

x(t,o,u) = mM ^didr^did}f^dg_ 

du dt* du dh* du du 

_ d£{^^ dg{t^ 



where we used the fact that 



dh* 8 7 
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Provided that t* is a smooth function, it can be expanded for h = around the 
real critical point T = and u = Uc, 

t*{t,0,u) = aifi + a2T + -- - , (2.20) 

with /i = M — Mc and r = T — T^. The absence of constant terms in the above 
expansion is due to Eq. f l2.15p because Aa{t*{tc, 0, Uc)) = fixes the temperature in 
the constrained system so that t*{tc,0,Uc) = 0. 
Then to first order 

+ (2.21) 
which inserted into Eq. f l2.18p with Eq. f l2.16p gives 

x(t, 0, u) = aiAo + aiA\t*\ + aiEiri^-" + "^^^^'^'""^ + ■ ■ • . (2.22) 

ou 

We can also expand the constraint, Eq. fl2.15p . about the real critical point u = 
and T = Tc 

X{t, 0, u) = X{tc, 0, Uc) + difi + d2T + -- - . (2.23) 
In addition, from Eq. (E^O]) 

^ = Lt*it,0,u)-—T+---. (2.24) 
ai ai 

Using the above equation, we can insert Eqs. fl2.22p and (12.231) into (12.121) to obtain 
the main result 

al{A\t*\ + 5|r|^-") = dit* + (aid2 - ^^102)^ . (2.25) 

If a < 0, the regular term dominates and oc |r|, resulting in the absence of 
Fisher renormalization, in which case the critical exponents of the transition remain 
unchanged. On the contrary, if a > 0, one obtains the central result 

|r| oc |r|^/(^-"). (2.26) 

That means that a deviation from the critical point of r in the real system is equival- 
ent to a deviation t* in the ideal system, these deviations being related by Eq. f l2.26p . 
Then the real system approaches the transition more slowly than the ideal one. The 
internal energy of the real system for /i = is 

df(t,0,u) ,dt*(t,0,u) dg(t,0,u) , 

e(t,0,n) = '^g^' ^ =eo(t*,0) ' ^ + ^ (2.27) 

= (Ao + A|r|+5|tt-")^^^%^ + ^%^ (2.28) 

at at 

= (Ao + A|r|Va-") + B\r\fJ^^^^ + , (2.29) 

at at 
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where the only possible singular contributions to the specific heat stem from the 
terms within parentheses. The specific heat of the constrained system is then 

c(t,0,«) = ^^^|^~|rr/(-"), (2.30) 

and therefore the singular behaviour of the specific heat has been replaced for a 
cusp-like one, i.e., it remains finite due to the presence of the constrained hidden 
variable. Then the Fisher renormalization for the exponent of the specific heat is 

ax = . (2.31) 

We can obtain analogously the renormalization for other critical exponent as: 

a ~ iri'^- |r|^/(i-") , /3^ = ^— , (2.32) 

1 — a 

X ~ |r|-^~|rr^/(i-") , 7x = ^, (2.33) 

I — a 

^ ^ It*]-" M-'/^'^"^ , ux = -^. (2.34) 

1 — a 

The critical exponent r] is defined at just the critical point and therefore is invariant 

Vx = V- (2.35) 

If the standard power-law scaling behaviour of the main quantities is modified 
by multiplicative logarithmic correction, i.e, 

Co(t) ~ |tr"|log|t|r, (2.36) 
mo(t) ~ |t|^|log|t||'^ fort<0, (2.37) 
Xo{t) ~ |tr^|log|t|r, (2.38) 

Co{t) ~ |tniog|tir, (2.39) 

mo{h) ~ \h\^\og\h\f for t = , (2.40) 
goix,t) ~ a:-(^-2+,)^jQg^),)^|^_|_^ fort«l, (2.41) 

then Eq. ( I2.10p is naively changed to give: 

/o(t, 0) = A±+A±N+^2±iti'+0(iti')+s±|t|2-"i log |tir |i + o 



log 


log 




log 


ltl 





(2.42) 

And therefore it can be obtained the equivalent of Eq. fl2.25p for the logarithmic 

case 

al{A\t*\ + 5|r log |r 11°) = dit* + (airfa - (^1^2)^ • (2.43) 
The above equation produces the following results: 



^ This can also be obtained from the scahng relationship (2 — i])!/ = 7. If and 7 are renorm- 
alised, rj remains unchanged. 
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• If a < 0, or a = and a < 0, the regular term dominates and \t*\ oc |r|, 
leading to the absence of Fisher renormalization. 

• If a > 0, or a = and a > 0, one obtains the modified central result 

\t*\ oc |r|i/(i-")|log|r||-"/(i-") , (2.44) 

which produces the renormalization of the individual exponents of the logar- 
ithms: 

ax = -j^, (2.45) 
1 — a 

Px = (2.46) 
i — a 

IX = 7 + (2.47) 
1 — a 

z>x = z>+r^. (2.48) 
1 — a 

Again, no renormalization takes place for the exponents fjx = fj and 6x = S . 
2.2.2 The Microcanonical Ensemble 

The first step in the construction of the ensemble is an extension of the configuration 
space. We add N{= L^) real momenta, pj, to our original variables, CTj (named 
spins here) P^IMj . Note that this extended configuration, {crj,^,}, appears in many 
numerical schemes (consider, for instance. Hybrid Monte Carlo [SU] simulations in 
Lattice Gauge Theory). We shall work in the microcanonical ensemble for the {cxj, pj} 
system. 

Let U be the original spin Hamiltonian (i.e., Eq. f l2.82p in our case). Our total 
energy is E| 

N 2 

£ = ^^ + U {e = £/N , u = U/N) . (2.49) 

i=l 

The momenta contribution, 

N 2 

Nk^Y.^, (2.50) 

1=1 

is necessarily positive, and it is best thought of as a "kinetic" energy. In this mech- 
anical analogue, the original spin Hamiltonian U can be regarded as a "potential" 
energy. 



^Note that this microcanonical ensemble exactly matches the conditions in the original Fisher 
work [58] : the momenta are some hidden degrees of freedom in thermal equilibrium with the spins, 
and a global constraint is imposed. It is also curious to rederive the results in Sec. l2.2.2] considering 
r momenta per spin (in this work F — 1, while Lustig JM; always considered F — 3). If one takes 
the limit F ^ oo, at fixed N, the canonical probability is recovered for the spins. 
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The canonical partition function is (/3 = 1/T) 

ZnW = / U^p J: e-^^ = (I) e-^\ (2.51) 

where Y2{a } denotes summation over spin configurations. Hence, the {pi} play the 
role of a Gaussian thermostat. The {pi} are statistically uncorrelated with the spins. 
Since {Ky^-^"^ = 1/(2/3), one has le)f~^^ = {u)f~^^ + 1/(2/3) . 

Furthermore, given the statistical independence of n and u, the canonical probab- 
ility distribution function for e, P^\e), is merely the convolution of the distributions 
for K and u: 

I'OO 

Pf\e)= dKPf^^\K)Pf^'\e-K). (2.52) 

^0 

In particular, note that for spin systems on a finite lattice, P^^'*'"('u) is a sum of (order 

N) Dirac 5 functions. Now, since the canonical variance of k is l/(/3-\/2iV), roughly 
y/N discrete m- levels, with u ^ e — 1/(2/3), give the most significant contribution 
to Pl^\e). We see that the momenta's kinetic energy provides a natural smoothing 

of the comb-like P^^^'^{u). Once we have a conveniently smoothed P^\e)., we may 
proceed to the definition of the entropy. 

In a microcanonical setting, the crucial role is played by the entropy density, 
s(e, A^), given by 

exp[Ns{e,N)]= JJ dp^ ^ 5(Are - £) . (2.53) 

^=l 

Integrating out the {pi} using the Dirac delta function in f l2.53p we get 



JV 
2 



exp[iV.(e,iV)] = -^^J]a;(e,n,iV), (2.54) 



JV-2 



u{e,u,N) = {e-u) — e{e-u), (2.55) 

where P is the gamma function and the step function, 6{e—u), enforces e > u. Equa- 
tion (12.541) suggests defining the microcanonical average at fixed e of any function 
of e and the spins, 0(e, {cTj}), as [51] 



y^/„ \ 0(e, |(Tj|) u;(e, u, N) 
{0)e ^ V ' ^ (2.56) 



We use Eq. fl2.54p to compute ds/de 

ds{e,N) 



(/3(e;{a,}))e, (2.57) 



de 
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Bearing in mind the crucial role of the generating functional in Field Theory (see 
e.g. [Z]), we extend the definition fl2.53p by considering a linear coupling between 
the spins and a site dependent source field hf. 



poo N 

exp[Ns{e, {hi}, N)] = JJ dpi ^ e^» ^'""^diNe - £) 



(2.59) 



where S = Ne is still given by Eq. (12. 49 p . without including the source term. In 
this way, the microcanonical spin correlation functions follow from derivatives of 
s(e,{/i,},iV): 



d[N. 



dhk 
d'^[N s] 



dhkdhi 



e,{h,},N 



{o'k)e,{h,} , 



{o-kCri)e,{h,} - {o-k)e,{h,} {cri)e,{h,} 



(2.60) 



In particular, if the source term is uniform hi = hwe observe that the microcanonical 
susceptibility is given by standard fluctuation-dissipation relations, see Ref. [7] and 
Eq. (1239]) below. 



Ensemble equivalence 

Equation (I2.53p ensures that the canonical probability density function for e is 

^^(e) = ^^exp[N{s{e,N)-Pe)] , (2.61) 

hence, Eq. iKET} . 

\0gPf\c2) - \ogPf\c,) = N r de ((/3)e - /?) , (2.62) 

where log means natural logarithm everywhere in this work. 

The relation between the canonical and the microcanonical spin-values is given 

by 

/oo 
dc{0)ePf\c). (2.63) 
-00 

Now, Eqs. fl2.6ip and fl2.63p imply that the canonical mean value will be dominated 
by a saddle-point at e^^, 

(/3)e|p, = /3 , (2.64) 

which can be read as yet another expression of the Second Law of Thermodynamics, 
Tds = de. 

The condition of thermodynamic stability (namely that (/3)e be a monotonically 
decreasing function of e) ensures that the saddle point is unique and that e^^ is 
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a maximum of Pp{e). Under the thermodynamic stabihty condition and if, in the 
large L hmit, 

d(/3)e 



de 



<0, (2.65) 



the saddle point approximation becomes exact: 



and we have ensemble equivalence: 



It follows that the microcanonical estimator 

CjL,e)= ^ , (2.68) 

evaluated at g^lL^ ^ will tend in the large-L limit to minus the canonical specific 
heat. Thus, if the critical exponent a is positive, Eq. f l2.65p will fail precisely at Cc. 
Hence, Eq. fl2.67p can be expected to hold for all e but Cc (or for all (3 but /3c). 

Double-peaked histogram 

The situation can be slightly more complicated if Pj3X^) presents two local maxima, 
remindful of phase coexistence. This is actually the case for one of our models 
- the D = 2, four-state Potts model [81]. From Eq. f l2.62p it is clear that the 
solution to the saddle-point equation fl2.64p will no longer be unique. We borrow 
the following definitions from the analysis of first-order phase transitions (where 
true phase coexistence takes place) [13]: 

• The rightmost root of Eq. fl2.64p . ^ , is a local maximum of P^^'' correspond- 
ing to the "disordered phase" . 

• The leftmost root of Eq. (12.641) . ^ , is a local maximum of Pj^'' corresponding 
to the "ordered phase" . 

• The second rightmost root of Eq. (12.641) . ^ , is a local minimum of P^^''. 

Maxwell construction yields the finite-system critical point, /3c,l, see Fig. 12.91 and 
Appendix [El 

0= //''"^de ((/3)e-/3c,L) , (2.69) 
and the finite-system estimator of the "surface tension" 



2L^-i 



/"/''"^de ((/3)e - /3c,l) . (2.70) 



Of course, in the large-L limit and for a continuous transition, E — )■ 0, — )■ /3, 
and e^^^^ ^, e^^^^ ^ — )• Cc , as we will see. 
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2.2.3 Our Microcanonical Finite-Size Scaling Ansatz 

Usually, the Microcanonical FSSA takes an entropy density scaling form [6TI463] . 
In close analogy with the canonical case, one assumes that s{e, {hg}, N) can be 
divided into a regular part and a singular term Ssing(e, {hg},N). The regular part 
is assumed to converge for large L (recall that = L^) to a smooth function of 
its arguments. Hence, all critical behaviour comes from Ssmg{e, {h^}, N). Note as 
well that we write {h^}, instead of {/ij}, to emphasise the spatial dependence of the 
sources (supposedly very mild [7]). Hence, 

s,,,,,{e,{h,},N) = L-''g(^L^{e-e,),{Ly^h,}') . (2.71) 

Here, g is a very smooth function of its arguments, while i/h = 1+^^^ is the canonical 
exponent, see e.g. [7], which does not get Fisher-renormalised. Corrections to FSS 
due to irrelevant scaling fields, have not played a major role in several previous 
analysis [6TH63] (in [63] only analytical scaling corrections were considered), but 
will be important for our precision tests. Leading order corrections were, however, 
explicitly considered in Ref. 

We will propose here alternative forms of the ansatz f l2.7ip . more suitable for a 
numerical work where neither Cc nor the critical exponents are known beforehand. 

Our first building block is the infinite-system microcanonical correlation length, 
^oo,e ■ Indeed, ensemble equivalence implies that, in an infinite system, the long- 
distance behaviour of the microcanonical spin-spin propagator G{f; e) = {<7xCrx+r)e — 
(o"x)e(c"x+r)e behaves for large f^as in the canonical ensemble (close to a critical point 
^oo,e is large, so that rotational invariance is recovered in our lattice systems): 



A 



G(r; e) = — ^e-^/«- , (2.72) 



where A is a constant. In particular, note that ensemble-equivalence implies that the 
anomalous dimension rj does not get Fisher-renormalised. We expect ^oo,e = 
if the correspondence between e and T are fixed through e = (e)'^^^'^^^^ . 

The basic assumption underlying the FSSA is that the approach to the L — )■ oo 
limit is governed by the dimensionless ratio L/^oo,e- Hence, our first form of the 
microcanonical FSSA for the observable O whose critical behaviour was discussed 
in referring to Eq. ( 12. 2 p is 

{0)L,e = L-^fo{L/U,e) + ■■■ . (2.73) 

In the above, the ellipsis stands for scaling-corrections, while the function fo is 
expected to be very smooth (i.e., differentiable to a large degree or even analytical). 
A second form of the microcanonical FSSA is obtained by substituting the scaling 
behaviour ^oo,e oc |e — ed"^™: 

{0)L,e = L-^fo {L'/'-ie - ee)) + ■ • • . (2.74) 
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Again, fo is expected to be an extremely smooth function of its argument @. In par- 
ticular, this is the form of the ansatz that follows from Eq. f l2.7ip by differentiating 
with respect to e or from the source terms. 

However, the most useful form of the microcanonical FSSA is obtained by apply- 
ing Eq. fl2.73p to the finite-lattice correlation length ^^..e, obtained in a standard way 
(see Ref. [7]) from the finite-lattice microcanonical propagator. We expect ^L,e/ L to 
be a smooth, one-to-one function of L/^oo,e; that can be inverted to yield L/C,oo,e as 
a function of C,L,e/L. Hence, our preferred form of the FSSA is 



{0)L,e = L — 



o I —r- I + L '^Go 



L \ L 



(2.75) 



Here, Fq and Go are smooth functions of their arguments and u is the first universal 
scaling corrections exponent. 

It is important to note that exponent u does not get Fisher-renormalised. Indeed, 
let us consider an observable O with critical exponent xq at a temperature T such 
e = {e)f^'^]T^. Now, ensemble equivalence tells us that Of^^^'^^ = OL=oo,e and that 
^f=^',T^ = ^L=co,e- Ehminating T in favour of see e.g. [7], we have 

0—1 ^ ^^o/. + Bo^l^^^, + ■••], (2.76) 

where Aq and Bq are scaling amplitudes. It follows that = cu, and that xq/i^ = 

XO,ni/ ^ni- 

The Quotient Method 

Once we have Eq. 02.751) . it is straightforward to generalise the quotient method |70] . 
In Appendix [B] we also describe how it should be modified in the presence of (mul- 
tiplicative) logarithmic corrections to scaling. 

Let us compare data obtained at the same value of e for a pair of lattices Li = L 
and L2 = sL with s > 1. We expect that a single e^Li^Lj exists such that the 
correlation-length in units of the lattice size coincides for both systems: 

L ~ sL ■ ^ 

Hence, if we compare now in the two lattices the observable O in (12. 75 p . precisely 
at ec,L,sL, we have 

^,'^"'"'"'" =3^ [I + Ao,sL- + ■■■], (2.78) 

where Ao,s is a non-universal scaling amplitude. One considers this equation for 
fixed s (typically s = 2), and uses it to extrapolate to L = 00 the L-dependent 
estimate of the critical exponents ratio xo,m/^m- At the purely numerical level, it 



''Note that the microcanonical weight (|2.55p is not analytical at each energy level of the spin 
Hamiltonian. 
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needs to be noted as well that there are strong statistical correlations between the 
quotients in fl2.77p and in ( 12.78 p . that reduces the statistical errors in the estimate 
of critical exponents. These errors can be computed via a jack-knife method, see 
e.g. [7]. 

In this chapter, we shall compute the critical exponents from the following op- 
erators (x is the susceptibility, while ^ is the correlation length, see Sec. 12.31 for 
definitions): 

X ^ xo = z/m(2 - 77) , (2.79) 

de^ Xo = iym + l. (2.80) 

The L dependence of ec,L,s follows from Eq. fl2.74p as applied to ^l/L for the two 
lattice sizes L and sL [7|l55]: 

ec,L,. = ec + — + ■ ■ ■ , (2.81) 

where B is again a non-universal scaling amplitude. In particular, if one works at 
fixed s, CcL sL tends to Cc for large 



2.3 The Model 

We will define here the model and observables of a generic Z)- dimensional Q-state 
Potts model. The numerical study was done for two instances of this model: the 
three-dimensional Ising (Q = 2) model, and the two-dimensional Q = 4 Potts model. 

We place the spins ai = 1, . . . ,Q at the nodes of a hypercubic D-dimensional 
lattice with linear size L and periodic boundary conditions. 

The Hamiltonian is 

W = -5^5.,., , (2.82) 

(id) 

where {i,j) denotes first nearest neighbours and 6ij is the Kronecker delta. For a 
given spin, a, we define the normalised Q-vector s, whose g-th component is 




A Q components order parameter for the ferromagnetic transition is 

A^ = -^J]sl, (2.84) 



where i runs over all the lattice sites. We will now consider microcanonical averages. 
The spatial correlation function is 



C{r'-r) = (s(r)-s(r') 

g / 1\ (2-85) 



^Note that, Eq. (|2.74p tells us that, if the energy histogram is double-peaked, see Sec. 12.2.21 the 
histogram maxima will tend to Cc only as L"^/''" . 
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Our definition for tlie correlation lengtli at a given internal energy density e, is 
computed from the Fourier transform of C 

C{k) = Y,C{r)^'^-\ (2.86) 

r 

at zero and minimal (||fcmin|| = 27r/L) momentum [7l[i 



2 sm(7r/L) 

Note that C can be easily computed in terms of the Fourier transform of the spin 
field, s(fc), as 

(7(fc)=L^(s(fc)-S(-fc)>^, (2.88) 
and that the microcanonical magnetic susceptibility is 

X = L^{M^)e = C{^). (2.89) 

For the specific case of the Ising model, the traditional definitions, using = ±1 
(recall that Si = ±l/-\/2), are related with those of the general model through: 



Rising ^ S,Sj =2U- 3L^ 

<i,j> 



Rising _ ^/2, (2.90) 

Notice that for D = 2 this model undergoes a phase transition at Pc = log(l + 
/Q) which is of second order for Q < 4 and first order for Q > 4 |82] . 



2.4 Numerical Results 



2.4.1 Methods 

We have simulated systems of several sizes in a suitable range of energies (see 
Table 12. ip . To update the spins we used a Swendsen-Wang (SW) version of the 
microcanonical cluster method [13]. This algorithm depends on a tunable para- 
meter, K, which should be as close as possible to {f3)e in order to maximise the 
acceptance of the SW attempt (SWA). This requires a start-up using a much slower 
Metropolis algorithm for the determination of k. In practice, we performed cycles 
consisting of 2 x 10^ Metropolis steps, a k refresh, 2 x 10^ SWA, and a further k 
refresh. We require an acceptance exceeding 60% to finish these pre-thermalization 
cycles fixing k, for the following main simulation, where only the cluster method is 
used. 

In both cases studied, we observed a very small autocorrelation time for all 
energy values at every lattice size. In the largest lattice for the four-state Potts 
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model we also considered different starting configurations: hot, cold, and mixed 
(strips). Although the autocorrelation time is much smaller, for safety we decided 
to discard the first 10% of the Monte Carlo history, using the last 90% for taking 
measurements. 



Model 


L 


A^m(xl06) 


A^c 


Energy range 


Q = 2, D = 3 


8 


20 


42 


[-0.8,-0.9] 




12 


20 


42 


[-0.8,-0.9] 




16 


20 


49 


[-0.8,-0.9] 




24 


20 


25 


[-0.845, -0.875] 




32 


20 


16 


[-0.87, -0.860625] 




48 


20 


10 


[-0.87, -0.860625] 




64 


5 


10 


[-0.870625, -0.865] 




96 


5 


10 


[-0.870625, -0.865] 




128 


5 


7 


[-0.869375, -0.865625] 


g = 4, D = 2 


32 


1024 


61 


[-1.2,-0.9] 




64 


128 


61 


[-1.2,-0.9] 




128 


32 


41 


[-1.08,-0.98] 




256 


32 


24 


[-1.08,-1.005] 




512 


25.6 


32 


[-1.07,-1.01] 




1024 


6.4 


30 


[-1.06,-1.02] 



Table 2.1: Simulation details for the two models considered. For each lattice size 
L we show the number of measurements A/'m at each energy and the total number 
of simulated energies uniformly distributed over the displayed energy range Ae. For 
the Q = 4, D = 2 model, the values of reported were reached only at specific 
energies near the peaks of the Maxwell construction, where additional energy values 
were simulated. 



2.4.2 D = 3 Pure Ising Model 

In Fig. 12.11 (upper panel) we show a scaling plot of the correlation length (in lattice 
size units) against (e — CcjL^^'^"'- For the susceptibility we plot x ~ L'^~^ (lower 
panel). If the data followed the expected asymptotic critical behaviour with mi- 
crocanonical critical exponents they should collapse into a single curve. In Fig. 12.11 
we have used the canonical critical quantities from Refs. [HSIIBI] transformed to the 
microcanonical counterparts using Eq. (12. 2p . From the plot it is clear that important 
scaling corrections exist in both cases for the smallest lattices, although they are 
mainly eliminated in the largest systems. 

To obtain the microcanonical critical exponents we used the quotient method, 
see Sec. 12.2.31 The clear crossing points of the correlation length for different lattice 
sizes can be seen in Fig. 12.21 The determination of the different quantities at the 
crossings, and the position of the crossing itself, requires one to interpolate the 
data between consecutive simulated energies. We found that the method of choice. 
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-0.5 





(e - e„) L 



0.5 



Figure 2.1: Scaling plot of the correlation length (in lattice size units) and the 
scaled susceptibility for the three-dimensional Ising model. We used the critical 
values, Cc = —0.867433 and z/m = 0.7077. Notice the strong scaling corrections for 
small systems as well as the data collapse for the largest ones. 



given the high number of energy values available, is to fit, using the least squares 
method, a selected number of points near the crossing to a polynomial of appropriate 
degree. Straight lines do not provide good enough fits. However, second and third 
order polynomials give compatible results. In practice, we fitted a second-order 
polynomial using the nine points nearest to the crossing, also comparing the results 
with those using the seven nearest points that turn out to be fully compatible. For 
error determination we always used a jack-knife procedure, see Appendix O 



L 


ec,L,2L 








8 


-0. 


861831(12) 


0.44922(3) 


0.8033(42) 


0.0564(2) 


12 


-0. 


865010(10) 


0.46106(5) 


0.7968(31) 


0.0492(4) 


16 


-0. 


866020(6) 


0.46710(5) 


0.7717(22) 


0.0469(4) 


24 


-0. 


866767(3) 


0.47411(4) 


0.7665(11) 


0.0437(3) 


32 


-0. 


867034(4) 


0.47813(6) 


0.7594(13) 


0.0425(5) 


48 


-0. 


867228(2) 


0.48278(5) 


0.7492(5) 


0.0412(3) 


64 


-0. 


867302(2) 


0.48555(11) 


0.7457(16) 


0.0397(8) 



Table 2.2: Lattice size dependent estimates of critical quantities for the microcanon- 
ical D = 3 Ising model. The displayed quantities are: crossing points ec,L,2L for the 
correlation length in units of the lattice size, ^/L itself at those crossing points, and 
the estimates of the correlation length exponent u^^ and the anomalous dimension 
rj^. All quantities were obtained using parabolic interpolations. 
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The numerical estimates for Cc, ^L,ec/L and the critical exponents u^a and rj^, 
obtained using the quotient method for lattice pairs (L, 2L) are given in Table 12.21 
Our small statistical errors allow one to detect a tiny L evolution. An extrapolation 
to infinite volume is clearly needed. 

Before continuing, let us recall our expectations as obtained by applying Fisher 
renormalization to the most accurate determination of canonical critical exponents 
known to us [ura = z^/(l — «) = i^/ {Du — 1)]: 

= 0.7077(5) (from u = 0.6301(4) [8]) , 
rj^ = r] = 0.03639(15) [85] , 
CO = 0.84(4) [g . 

Besides, although non-universal, let us take Cc = —0.867433(12) [^. 

The results obtained from an extrapolation using only leading order 
rections were: 

• Cc = -0.867397(6), cu + = 1.918(26) 
(we obtained a good fit for L > L^^^ = 12, with xVd-o.f. = 0.39/3, C.L.=94%, 
where "d.o.f." stands for degrees of freedom and "CL." for confidence level^. 

• = 0.5003(12), CO = 0.581(27) 

(Li^in = 12, xVd.o.f. = 0.12/3, C.L.=99%). 

• Ura = 0.714(28), CO = 0.53(30) 

(L„,in = 8, xVd-o.f. = 3.16/4, C.L.=53%). 

• r]m = 0.0391(15), CO = 1.21(24) 

(L^in = 8, xVd.o.f. = 0.96/4, C.L.=92%). 

The main conclusions that we draw from these fits are: (i) the exponents are com- 
patible with our expectations from Fisher renormalization, (ii) sub-leading scaling 
corrections are important given the tendency of the fits to produce a too low estimate 
for 00 (see below), and (iii) the estimates from canonical exponents (themselves ob- 
tained by applying the high-temperature expansion to improved Hamiltonians [H|I55] ) 
are more accurate than our direct computation in the microcanonical ensemble. 

We can, instead, take an opposite point of view. If we take the central values in 
Eqs. fl2.91[ [2^921 12.931) as if they were exact, we can obtain quite detailed information 
on the amplitudes for scaling corrections: 

• We find an excellent fit to Um{L, 2L) = z/^ + AiL'"^ + AsL^^"^, for Lmm = 16: 
xVd.o.f. = 1.53/3, C.L.=68%, with = 1.38(7) and A2 = -7.6(1.1). This 
confirms our suspected strong sub-leading corrections. Indeed, according to 
these amplitudes Ai and A2, only for L ^ 130 the contribution of the (sub- 
leading) quadratic term becomes 10% of that of the leading one. 

*^For the 3D Ising model at criticality, wi'''"s = -0.990627(24) [83], and ^J"'"s = 0.2216546(2) [84], 
we obtain for our Potts representation of the Ising model Cc = (wj."™^ — D)/2 + l/{Af3l^™^). 

^The confidence level is the probability that would be larger than the observed value, sup- 
posing that the statistical model is correct. As a rule, we consider a fit not good-enough whenever 
C.L.< 10%. 



(2.91) 
(2.92) 
(2.93) 

scaling cor- 
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Figure 2.2: Crossing points of the correlation length in lattice size units for the 
three-dimensional Ising model. The error bars are in every case smaller than the 
point sizes. The values of the different quantities at the crossing as well as the 
critical exponents are given in Table 12.21 



• In the case of r]ra{L, 2L) = r]^ + BiL^'^ + 321-^"^ , for Lmin = 8: xV^-O-f. = 
2.4/5, C.L.=79%, we have = 0.101(10) and B2 = 0.07(7). Sub-leading 
scaling corrections are so small that, within our errors, it is not clear whether 
or not B2 = 0. 

The quite strong scaling corrections found for may cast some doubt on the 
extrapolation for C,L,ec/L, the only quantity that we cannot double check with a 
canonical computation. To control this, we proceed to a fit including terms linear 
and quadratic in with u = 0.84(4). We get 

%^ = 0.4952(5)(7), 

with Lmin = 12, x^/d.o.f. = 2.17/3, C.L.=54%. Here, the second error is due to 
the quite small uncertainty in u. It is remarkable that the contribution to the error 
stemming from the error in u is larger than the purely statistical one. 

The canonical specific-heat 

Previous numerical studies of microcanonical FSS focussed on the specific 

heat. Although we show all across this paper that a complete microcanonical FSS 
analysis can be based only on the spin propagator, the specific heat can be certainly 
studied within the present formalism. 
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As discussed in Sect. 12.2.21 (see also [13]), the canonical specific heat can be 
estimated from the microcanonical estimator C^{L,e) defined in Eq. fl2.68p . The 
expected FSS behaviour for C^iL, ec,L,2L) is 

C^{L, ee,L,2L) = i^'/^Ao + AiL-'^ + ■■■] + . (2.94) 

Here, Aq and Ai are scaling amplitudes, while i? is a constant background usually 
termed analytical correction to scaling, stemming from the non-singular part of the 
free-energy [7]. It is usually disregarded as it plays the role of a subleading scahng- 
correction term. Yet, a peculiarity of the D = 3 Ising model is that B is anomalously 
large (see e.g. [BB]) and needs to be considered. 

In Fig. 12. 3[ we reproduce the analysis of Bruce and Wilding [63], where the 
amplitude Ai in Eq. fl2.94p was fixed to zero by hand. In this way, if we consider 
the range of lattice sizes 8 < L < 64 (in [HB] only L < 32 was considered), we 
obtain B = —35.01(11) but with an untenable x^/d.o.f = 227/5. Our value of B is, 
nevertheless, quite close to the result B = —34.4(4) reported in [63] (unfortunately, 
these authors provided no information on fit-quality). 

Once the arbitrary constraint = is removed, we do obtain an acceptable 
fit, x^/d.o.f = 0.68/4. Perhaps unsurprisingly, the estimate of B is largely changed, 
once a nonvanishing Ai is allowed: B = —24.4(7). 




Figure 2.3: Microcanonical estimate of the specific heat, Cm{L,e), at ec,L,2L for the 
D = 3 Ising model, as a function of the system size. The numerical estimates of 
exponents a/u and u were taken from Ref. p]. The error bars are in every case 
smaller than the point sizes. The solid line is a fit to Eq. fl2.94p (fitting parameters: 
Aq, Ai and B), the dashed one is obtained by constraining the fit to Ai = 0. 
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2.4.3 D = 2, Q = 4 Pure Potts Model 

The Q = 4, D = 2 Potts model involves two peculiarities that will be explored here. 
First, it suffers from quite strong logarithmic scaling corrections. And second, it 
displays pseudo-metastability |81], an ideal playground for a microcanonical study. 

The study of the FSS for the Q = 4, D = 2 Potts model based on the 
analysis of the Renormalization Group (RG) equations ^H], reveals the presence of 
multiplicative logarithmic scaling corrections. This is one of the possible forms that 
scaling corrections can take in the limit w — )■ 0, and is a major nuisance for numerical 
studies. A very detailed theoretical input is mandatory to safely perform the data 
analysis. We shall make here an educated guess for the microcanonical form of the 
scaling corrections, based purely on ensemble equivalence and the canonical results. 

From ensemble equivalence we expect 

e-ee~C(L,/3e)/\/3L , (2.95) 

where C{L, /3c) is the finite-lattice canonical specific heat at /5c, and A/3 = /3c^^ — /3c 
is the inverse-temperature distance to the critical point of any L-dependent feature 
(such as the temperature maximum of the specific heat, etc.). We borrow from 
Ref. |13] the leading FSS behaviour for these quantities: 

Thus, we have: 

e(L) - e,(oo) ~ L-'/^logL)-^/'' . (2.97) 

This result can be derived as well by considering only the leading terms of the 
first derivative of the singular part of free energy with respect to the thermal field, 

(oc /3 - /3c) 



^^^'"^^'^^'^''^^ - Id^\<P\'^\- log 101)-^ + D^\<pm- log 101)-^^ . (2.98) 
The above equation describes the energy of the system, and its leading term is 

e-ec~-D±f^, 2.99 
3 log 101 

but 

(l)^C'^L-^/^(logL)^/\ (2.100) 

so it is direct to obtain again Eq. fl2.97p . Hence, we are compelled to recast Eq. (12.741) 

as 

{0)L,e = L-^fo {L'/'i\ogLf\e - Cc)) +■■■. (2.101) 

Furthermore, from the canonical analysis |13], we expect multiplicative logarithmic 
corrections to the susceptibility (that do not get Fisher renormalised) . Furthermore, 
the ellipsis in fl2.10ip stands for corrections of order log log L / log L and 1 / log L [13] . 

We first address in the next subsection the direct verification of Eq. fl2.10ip using 
the quotient method. We then consider the pseudo-metastability features. 
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Scaling plots and critical exponents 




(e - e^) L (log L) 



Figure 2.4: Graphical demonstration of Eq. fl2.10ip as applied to the microcanonical 
D = 2, Q = 4 Potts model: both the correlation length in units of the lattice size 
(top) and the scaled susceptibility, x Eq. fl2.102p (bottom), are functions of the 
scaling variable (e — ec)I/^/^(logL)^/^. 



We start with a graphical demonstration of Eq. ( 12.1011) : ^/L as a function of 
(e — Cc) L^^"^ {log L)^^'^ should collapse onto a single curve (the deviation will be larger 
for small L values due to neglected scaling corrections of order log log L / log L and 
1/logL) 1^. Similar behaviour is expected for the scaled susceptibility [i5] : 

Note that ^/L does not need an additional logarithmic factor. These expectations 
are confirmed in Fig. 12. 4[ especially for the larger system sizes (that are subject to 
smaller scaling corrections). 

We can check directly the importance of the multiplicative logarithmic correc- 
tions for the susceptibility by comparing x and x as a function of ^/L, see Fig. 12.51 
The improved scaling of x is apparent. We observe as well that the largest cor- 
rections to scaling are found at and below the critical point (around ^/L ^ 1.0). 

The scaling proposed for the susceptibility in Ref. [13] can also be checked from 
our values at ecL2L- Considering (our data is fully supportive of this 

point) we can plot log(x/L'^/^) versus log log L. We obtain a linear fit for the data 
with L > 64 with a slope —0.132(3) (x^/n.d.f. = 7.5/1), see the dashed line in 

^ We obtain the exact Cc in the thermodynamic hmit from (3c ~ log(l + \/Q) [57], and Uc = 
-(1 + [52] by applying = + l/(2/3e). 
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Figure 2.5: Comparison of the scaling for the naively scaled susceptibility xL~'^^^ 
(top) and for x, defined in Eq. f l2.102p . (bottom), as a function of the correlation 
length in units of the lattice size, for the microcanonical D = 2, Q = 4 Potts model. 



Fig. I2.4.3[ which can be compared with the expected value —1/8 [13]. The large 
value of x^/d.o.f. can be ascribed to the presence of higher order correction terms. 
In fact the whole scaling behaviour for the susceptibility is pI3] 

X^L">gL)-'»(l + A!^ + B^ + ...). (2.103) 

and we can use this form for a least-square fit. Fixing both the leading and the 
logarithmic exponents we estimate A = 0.80(7) and B = —0.48(3) using all the 
lattice sizes with x^/d.o.f. = 2.9/2, see the solid line in Fig. 12.4.31 Therefore our 
data set is fully supportive of the behaviour proposed in Ref. |13], including the 
subleading additive logarithmic corrections. 

We now proceed to the numerical computation of critical exponents. We shall 
use the quotient method, modified as described in Appendix [Bl From Fig. 12. 7^ we 
can see that the crossing points can be well obtained using parabolic interpolations 
of the nine points around the estimated crossing energies, as done in Sec. 12.4.21 
We checked that the results do not depend on the interpolating polynomial degree 
by comparing with interpolations using cubic curves. We also compared with the 
results obtained using only seven points around the crossing obtaining again full 
agreement. 

The critical exponents obtained are listed in Table 12.31 They may be compared 
with the exact ones [82] (z/ = 2/3, a = 2/3 and rj = 1/4): 

u^ = 2 ; r^ = r]^ = ^. (2.104) 
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Figure 2.6: Logarithmic scaling behaviour of the susceptibility at the critical point. 
The error bars are in every case smaller than the point sizes. The dashed line do 
not include the subleading additive terms of Eq. fl2.103p while the solid line do. 



Comparing with our computed exponents, we obtain an acceptable agreement. In 
the case of the microcanonical u exponent, i^m, after adding the correction for the 
quotient method in the presence of logarithms, the agreement is fairly good. We can 
see a clear trend towards the exact value for all the lattice sizes except the largest 
(2.5 standard deviations away), which is probably due to a bad estimate of the huge 
temperature derivatives of the correlation length. In the case of the microcanonical 
T) exponent, r/m, which must be the same as the canonical one, we can see clearly 
the tendency to the analytical value r/m = 0.25. We must stress the importance of 
adding the corrections described in Appendix [B] to the quotient method. 



Critical point, latent heat, and surface tension 

It has been known for quite a long time that the D = 2, Q = A Potts model on finite 
lattices shows features typical of first-order phase transitions [HI]. For instance, see 
Fig. \2.S\ the probability distribution function for the internal energy, Pi3{e), displays 
two peaks at energies Cd (the coexisting disordered phase) and Co (the energy of the 
ordered phase) separated by a minimum at e*. Of course, since the transition is of 
second order, Cc is the common large L limit of Cd, Co and e*. 

We discussed in Sec. 12.2.21 how the Maxwell construction is used to estimate the 
canonical critical point /3c,l, as well as Cd, Cq and the associated surface tension. 
This procedure is outlined in Fig. 12.91 The numerical results are listed in Table 12. 4[ 
where we can see that Pc,l is a monotonically increasing function of L continuously 
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Figure 2.7: Correlation length in lattice size units for the D = 2, Q = A Potts model. 
The values of the different quantities at the crossings for lattices L and 2L, as well 
as the corresponding estimates for critical exponents, are given in Table 12.31 The 
inset is a magnification of the critical region. 



approaching the analytical value f3c = log{l + ^/Q) = 1.0986122 . . . [87]. A jack-knife 
method [7] was used to compute the error bars for all the quantities in Table 12.41 

To perform a first check of our data, we observe that Pc,l is a typical canonical 
estimator of the inverse critical temperature. As such, it is subject to standard ca- 
nonical FSS, where the main scaling corrections come from two additive logarithmic 
terms |43j: 

{hgLf^f log log L , 1 \ 

From our data in Table EH we obtain ai = —0.44(7), 02 = —1.15(72), and 03 = 
2.28(26), and a good fit (L^in = 128, xVd-o.f. = 0.28/1, C.L.=60%). 

As for the L dependence of Cd and Co, we try a fit that considers the expected 
scaling correction terms |43j : 

ec,o,L - = a,L-'/\\ogL)-'/' (l + a^^^f^ + a,-^] . (2.106) 

V logL logLy 

Our results for Cq are: = —2.03(20), = —1.65(27), and a^o = —2.08(41), 
with a fair fit quality (L^^in = 32, xVd-o.f. = 2/3, C.L.=57%). We obtain for 
^d- Old = 2.02(14), a2d = 0.93(37), and a^^ = —2.93(34) , with a fair fit as well 
(Lmin = 32, xVd-o.f. = 0.84/3, C.L.=84%). These two fits are shown in Fig. [2A01 

For the surface tension, one notes in Table 12.41 a non-monotonic behaviour. Fur- 
thermore, we lack any theoretical input with which to attempt a fit. We thus turn to 
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L 








< 






32 


-1.04659(5) 


0.8016(5) 


1.534(6) 


1.998(10) 


0.2663(9) 


0.2334(9) 


64 


-1.04633(2) 


0.7990(3) 


1.554(8) 


1.957(12) 


0.2638(6) 


0.2360(6) 


128 


-1.04579(1) 


0.7909(3) 


1.578(5) 


1.938(7) 


0.2639(5) 


0.2398(5) 


256 


-1.04548(2) 


0.7836(5) 


1.643(12) 


1.987(17) 


0.2615(11) 


0.2402(11) 


512 


-1.04519(2) 


0.7734(9) 


1.602(31) 


1.895(42) 


0.2617(21) 


0.2427(21) 



Table 2.3: Crossing points of the correlation length in lattice size units as a function 
of the energy for pairs of lattices (L, 2L). Using the original quotient method ^ we 
obtain the microcanonical critical exponents, listed in Columns 4 and 6, while the 
corrected ones (Columns 5 and 7) are labelled with primed symbols, see Appendix iBl 
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Figure 2.8: Canonical probability distribution function for the energy density, 
P^p^\e), as reconstructed from microcanonical simulations of the D = 2, Q = A 
Potts model for different system sizes. The L-dependent critical point /3c,l is com- 
puted using the Maxwell rule. Sec. 12.2.21 (note the equal height of the two peaks 
enforced by Maxwell construction). The system displays an apparent latent heat 
that becomes smaller with increasing L, and vanishes in the large L limit. 



a variant of the quotient method. Were S to follow pure power-law scaling, U oc L^, 
the exponent h would be obtained as: 




_ log[r(Li)/r(L2)] 
log(Li/L2) 



The effective exponent h obtained from our data is given in Table 12.51 We observe 
that it is clearly negative (as it should be since E vanishes for a second-order phase 
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Figure 2.9: Top: From the microcanonical mean values {P)e,L for the D = 2, Q = 4 
Potts model, we estimate the size dependent canonical inverse critical temperat- 
ure Pc,L (horizontal lines) for all the simulated lattice sizes, ranging from L = 32 
(lower) to L = 1024 (upper). We show as well the analytical prediction (uppermost 
horizontal line). Bottom-left: Example of Maxwell construction for our L = 32 
data. The e- integral of (/3)e,L — Pc,l from Co to Cd vanishes. Bottom-right: Zoom 
of upper panel showing only data for lattice sizes L = 256 (lower curve), L = 512 
(middle curve), and L = 1024 (upper curve). 



L 


/3c,L 


Co 




r X 10^ 


32 


1.0911070(20) 


-1.0175(4) 


-0.9760(2) 


0.47(2) 


64 


1.0957256(14) 


-1.0392(3) 


-0.9915(2) 


2.77(7) 


128 


1.0975150(10) 


-1.0463(3) 


-1.0062(5) 


4.10(15) 


256 


1.0981989(5) 


-1.0489(2) 


-1.0183(3) 


3.92(8) 


512 


1.0984570(3) 


-1.0490(1) 


-1.0266(2) 


3.28(11) 


1024 


1.0985539(3) 


-1.0483(3) 


-1.0325(1) 


2.09(17) 



Table 2.4: Using the Maxwell construction, we compute for the D = 2, Q = A 
Potts model the L-dependent estimates of the (inverse) critical temperature /3c,l, 
the energies of the coexisting ordered Cq, and disordered Cd phases, as well as the 
surface tension E. 



transition). An asymptotic estimate, however, seems to require the simulation of 
larger systems. 
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0.05 r-i/2 0.1 0.15 



Figure 2.10: System size dependent estimates of the energies of the "coexisting" 
ordered (co, blue squares) and disordered (cd, red circles) phases of the D = 2, 
Q = A Potts model, as a function of The lines are fits to the expected 

analytical behaviour Eq. f l2.106p . The horizontal line corresponds to the asymptotic 
value, Cc. 



(^1,^2) 




(32,64) 
(64,128) 
(128,256) 
(256,512) 
(512,1024) 


2.56(7) 
0.56(6) 
-0.065(60) 
-0.257(57) 
-0.650(127) 



Table 2.5: Effective exponent obtained using Eq. fl2.107p for the surface tension. 



We have just seen that, up to scaling corrections, e^^^ and ei^^ correspond 
to (different) L-independent values of the argument of the scaling function in 
Eq. f l2.10ip . Hence we expect that ^{ed)/L and ^{eo)/L, see Table I2l6| approach 
non- vanishing, different values in the large L limit. The FSS corrections are expected 
to be additive logarithms [43] 

^ = a + -\. (2.108) 



The results are: 



f(eo) , , 2.28(5) 

^V^ = 1.28 1 - , V , (2.109) 
L log-L 
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L 


aeo)/L 




x(eo) 


X(ed) 


^canonical ^ 


TTjcanonical 

A. 


32 


0.637(2) 


0.453(1) 


0.907(2) 


0.647(1) 


0.990(3) 


1.287(3) 


64 


0.732(3) 


0.396(1) 


1.025(3) 


0.545(2) 


0.995(2) 


1.310(2) 


128 


0.799(5) 


0.357(4) 


1.106(5) 


0.472(7) 


1.001(3) 


1.331(3) 


256 


0.866(6) 


0.335(3) 


1.182(6) 


0.429(5) 


1.001(5) 


1.343(5) 


512 


0.915(4) 


0.315(2) 


1.238(4) 


0.392(4) 


1.014(8) 


1.366(8) 


1024 


0.953(15) 


0.302(2) 


1.279(13) 


0.367(3)) 


0.997(21) 


1.353(22) 



Table 2.6: Correlation length in units of the lattice size and the RG invariant x 
defined in Eq. f l2.102p . for several L values, as computed in the microcanonical 
D = 2, Q = 4 Potts model. The values of the energy density correspond to the 
ordered (cq) and disordered (cd) phases. For comparison we also display in the last 
two columns the canonzca/ results at (3c obtained in Ref. [43j. 



(L^in = 32, xVd.o.f. = 4.2/3, C.L.=22%), and 

f(ed) , , 0.98(2) 

^V^ = 0.159 4 - - — V^, 2.110 
L log L 

(L^in = 32, xVd.o.f. = 3.3/3, C.L.=37%). 

A very similar analysis can be performed for the scaled susceptibility, Eq. fl2.89p . 
at Cd and Cq. In order to deal with the multiplicative logarithms of the susceptibility, 
we used x defined in Eq. fl2.102p . 

Fitting our data set to the logarithmic form 

__^^^lo^_ (2.111) 
logL 

obtained in Ref. [13], we obtain a good fit in the ordered phase energy, Co'. 

x(eo) = 2.41(5) - 4.00(15)i^f^ , (2.112) 

logL 

with Lmin = 128, x^/d.o.f. = 3.10/2, C.L.=21%. However, the extrapolation for the 
susceptibility defined in the disordered phase energy, Cd, is a nonsensical negative 
value. 

We can also fit the data to the logarithmic form also used in Ref. 



logL 

finding: 



X = A + -^, (2.113) 



logL 

(Linin = 32, xVd.o.f = 7.44/4, C.L. = 11%), and 



x(eo) = 1.643(5)-^^, (2.114) 



X(ed) = 0.094(7) + , (2.115) 

log 1j 
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(Lmm = 64, x^/d.o.f = 2.94/3, C.L.=37%) . For comparison, we recall that Ref. 
reports two different fits for x, depending on the logarithmic corrections they used: 

^anonical _ 1.673(33) - 1 .056 (98) ^"f ^ , (2.116) 

logL 

^anonical ^ 1.454(13) _ ^L^^^M . (2.117) 



2.5 Conclusions 

We have formulated the Finite Size Scaling Ansatz (FSSA) for microcanonical sys- 
tems in terms of quantities accessible in a finite lattice. This form allows to extend 
the phenomenological renormalization approach (the so-called quotient method) to 
the microcanonical framework. 

Our FSSA was subjected to strong numerical testing. We performed extens- 
ive microcanonical numerical simulations in two archetypal systems in Statistical 
Mechanics: the three-dimensional Ising model and the two-dimensional four-state 
Potts model. The two models present a power-law singularity in their canonical 
specific heat, implying non-trivial Fisher renormalization when passing to the mi- 
crocanonical ensemble. A microcanonical cluster method works for both models, 
hence allowing us to study very large system sizes [L = 128 in D = 3 and L = 1024 
in D = 2). 

In the case of the Ising model, we obtained precise determinations of the critical 
exponents that provide strong evidence for our extended microcanonical FSS ansatz. 

For the Potts model, very strong logarithmic corrections (both multiplicative 
and additive) plague our data. Fortunately, we have a relatively strong command 
over these corrections from canonical studies Our data can be fully rationalised 
using the scaling corrections suggested by the theoretical analysis ^3] . 



50 



Chapter 3 



Quenched Disorder Effect on a 
First-Order Phase Transition 

3.1 Introduction 

Although first-order phase transitions are by far the more frequent in nature, not 
much is known about the consequences of adding impurities to systems that in the 
pure case undergo this type of transition. This is due to the fact that there exist 
inherent difficulties for their study. 

One of the intrinsic problems in simulating first-order phase transitions is that 
in this case two or more phases coexist at the critical temperature. The system 
changes from the high temperature phase to the low temperature one by building 
an interface of size L, where L is the lattice size. The energy cost of such a mixed 
configuration is IJL^~^ (with E being the surface tension and D the spatial dimen- 
sion). Therefore, when doing simulations using the canonical ensemble (at fixed 
temperatures), the probability of reaching such mixed configurations is attenuated 
by a factor exp[— Z'L^~^], and as a result the natural time scale of the simulation 
grows with the system size L as exp[Z'L'^~^]. This huge obstacle to simulating large 
systems is called Exponential Critical Slowing Down (ECSD). 

Up to now, no solution for ECSD has been found in canonical simulations. This 
has motivated the popularity of simulations within the microcanonical ensemble (at 
fixed energy), see Sec. 12.2.21 Some simulation methods within this ensemble consider 
the canonical probability density function (pdf) of the energy as a constant within 
the energy interval e° < e < e'^ (e° and e*^ being the energy densities of the coexisting 
ordered and disordered phases respectively). This led to these methods being called 
flat-histogram methods [88H9T]. The canonical probability minimum in the energy 
gap (oc exp[— Z'L^~^]) is achieved by means of an iterative parameter optimisation. 
In fiat-histogram methods the system performs an energy random walk in the energy 
gap. The elementary step being of order (a single spin-fiip), one naively expects 
a tunnelling time from e° to e'^ of order spin-flips. But the (one-dimensional) 
energy random walk is not Markovian, and these methods still suffer ECSD [92] . In 
fact, for the standard benchmark (the Q = 10 Potts model j82j in D = 2), the barrier 
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of 10^ spins was reached in 1992 |SH], while the largest simulated system (to the best 
of our knowledge) had 4 x 10^ spins [89] . 

ECSD in flat-histogram simulations is probably understood [22]: on its way 
from e'^ to e°, the system undergoes several (four in Z) = 2) "transitions". First 
comes the condensation transition [921193] . at a distance of order L~^/*^^+^) from 
e'^, where a macroscopic droplet of the ordered phase is nucleated. Decreasing e, 
the droplet grows to the point that, for periodic boundary conditions, it reduces its 
surface energy by becoming a strip [M], see the figures in [T3| (in D = 3, the droplet 
becomes a cylinder, then a slab At lower e the strip becomes a droplet of 

disordered phase. Finally, at the condensation transition close to e° , we encounter 
the homogeneous ordered phase. 

In this work we will study a prototypical model of a strong first-order phase 
transition, the three-dimensional Potts model with Q > 3 states. There are numer- 
ous experimental systems which can be mapped by this model. For instance, the 
Q = 4 pure case in two dimensions describes the adsorption of N2 molecules on Kr 
in graphite layers |96j ; in three dimensions it describes the behaviour of FCC anti- 
ferromagnetic lattices (NdSb, NdAs, and CeAs, for example) with the magnetic field 
pointing in the (1,1,1) direction [97]. The site-diluted Q = 4 case in two dimen- 
sions models the effect of oxygen impurities on a sample of nickel where hydrogen 
molecules are adsorbed ^98j. In the dilute three-dimensional case we are not aware 
of any experimental realization. 

It is known fSlllB] that the pure three-dimensional Potts model undergoes a 
first-order phase transition in the pure case for Q > 3. On the contrary, it has been 
been found [H] that for strong dilution the system performs a second-order phase 
transition. A direct question is the following: what is the exact dilution that causes 
the order of the transition to change? What is more, are we absolutely sure that 
first-order phase transitions exist in the presence of dilution? This is still an import- 
ant open problem in Statistical Mechanics, and also one with implications in very 
technical fields such as highly correlated electron systems (e.g., high temperature 
superconductors or colossal magnetoresistance oxides) where phase coexistence and 
chemical disorder play crucial roles [99] . 

The question in the previous paragraph can be considered exactly solved in two 
dimensions [100] : even the most insignificant amount of impurities is enough to 
switch the phase transition from first-order to second-order (for the Universality 
Classes see [101] ). In D = 3 the most useful physical picture is provided by the 
Cardy-Jacobsen conjecture [101] : considering a ferromagnetic system undergoing a 
first-order phase transition for a pure sample, with T being the temperature and p 
the concentration of magnetic sites, a critical line, Tc{p), separates the ferromagnetic 
and the paramagnetic phases in the {T,p) plane. In D = 3 a critical concentration 
is expected to exist, 1 > Pc > 0, such that the phase transition is of first order for 
p > Pc and of second order for p < p^ (at p^ one has a tricritical point). When p 
approaches pc from above, the latent heat must vanish with the critical exponent 
of the magnetisation in the Random Field Ising Model (RFIM). Also the surface 
tension, E, vanishes at pc, while the correlation length ^(Tc(p)) diverges, with critical 
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exponents related to those of the The main objection to this argument is that 

the Cardy-Jacobsen conjecture relies on a mapping from the (large Q) disordered 
Potts model [82j onto the RFIM (two unsolved models in D = 3). As a result, if the 
D = 3 RFIM phase transition turned out to be of first order [104] . the conjecture 
would not be valid. 

The D = 3 problem has already been numerically studied [HHin] ; large regions 
of the critical line Tc{p) were found to be second order. Unfortunately, the study of 
the tricritical point as well as that of the first-order part of the critical line seemed 
beyond hope, mainly due to two factors. Firstly, an important difficulty arises from 
the long-tailed pdf 's encountered when comparing the specific heat or the magnetic 
susceptibility of different samples at Tc{p) Note that diverging- variance pdf's 
arise from the common practice of defining the quenched free energy at temperature 
T as the average of the sample's free energy at T [22], which is dominated by rare 
event^. Secondly, the other factor has been described above - the simulation of a 
sample of linear size L with previous methods is intrinsically difficult: the required 
simulation time grows exponentially with L^~^ [52] due to the ECSD. These two 
factors have limited previous work [ISllin] to L < 25. 

To overcome these two difficulties, on the one hand we propose two alternat- 
ive methods of performing the sample average, both of which reproduce the correct 
Thermodynamic Limit, avoiding the diverging-variance pdf's, and providing comple- 
mentary information, and on the other we exploit a novel microcanonic Monte Carlo 
method [13], which allows one to study the system entropy directly. This method, 
combined with a slightly modified typical cluster algorithm [TT|[T3]. permits accurate 
studies of systems with more than 10^ spins (when the previous methods can only 
handle 10'^). In our case the method will allow us to simulate systems of size up to 
L = 128 in the case Q = 4, D = 3, also making it possible to perform a Finite-Size 
Scaling (FSS) study of the elusive tricritical point as well as the associated critical 
behaviour. 

The highly accurate numerical study presented in this chapter has only been 
possible due to our capability of using different supercomputing facilities simultan- 
eously: 

• For the Q = A case: on the Mare-Nostrum machine of BSC (Barcelona Su- 
percomputing Centre) we used 160 000 computation hours (PowerPC 2.3 GHz 
processors); on the BIFI (Instituto de Biocomputacion y Fisica de Sistemas 
Complejos de Zaragoza) cluster we used 250 000 hours (Xeon Dual Core 3.40 

^ The expected exponents jS and v of the tricritical point |101| are: /3 = /3rfim and l/v ^ 
D- ^RFiM — /Jrfim/^'rfim Or (modified hyperscahng relation of the RFIM) v = j^rfim/ (2 — aRFiM^ 
Prfim)- The surface tension goes to zero with an exponent fi ~ {D ~ l)^ [102] , Taking the critical 
exponents of the Gaussian RFIM: /3rfim = 0.00(5), i^rfim — 1-1(2) and 0rfim = 1-53(1) [103] . the 
exponents for the tricritical point should be: ~ 1.5, /3 ~ and ^ ~ 3. 

^ Equilibrium phase-coexistence in a sample of N spins occurs for a temperature interval of 
width ^ iV~^ [105j . where the specific heat is C ~ A^. Yet the sample-averaged C scales at most 
as iV^/^ |106j because the sample dispersion of the critical temperatures leads to the critical region 
having a width ^ N^^/"^ around Tc{p). For any fixed temperature within the critical region, only 
a fraction ^ N~^^^ of the samples displays C ^ N. 
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GHz processors); and on computers (mostly Pentium 2.6 GHz) located in the 
UEX (Universidad de Extremadura) and UCM (Universidad Complutense de 
Madrid) we used 65 000 and 160 000 hours respectively. As a result we estim- 
ate that the computational resources used for this part are equivalent to 60 
years of a single last generation (Pentium 2.5 GHz) processor. 

• For the Q = 8 case: we used mostly the IBERCIVIS infrastructure, see Ap- 
pendix [Gj from which we obtained the huge number of approximately 300 years 
of a single last generation (Pentium 2.5 GHz) processor. In addition we made 
extensive use of the BIFI cluster, from which we obtained around 40 years of 
equivalent simulation time. We also used local resources in Badajoz but they 
can be disregarded compared to the aforementioned enormous numbers. 

3.2 Analytical Framework 

In this section we briefly review the main analytical results on first-order phase 
transitions with disorder. They have been taken from Ref. |100] . where it was 
demonstrated that for D < 2 even the smallest amount of impurities (whether in 
the bonds or in the fields) destroys the discontinuities of the first derivatives of 
the free energy making the transition continuous (of second order type), and from 
Ref. |101] where, after relating the dilute Potts model with the RFIM, it was found 
that for D > 2 there must exist a region in the phase diagram where the transition 
continues to be of the first order type even in the presence of disorder. This region 
will end up in a tricritical point. 

3.2.1 Aizenman-Wehr Theorem 

In Ref. |100] . it was demonstrated that for D <2 the presence of quenched random 
fluctuations in the structural parameters (external field h, temperature T, ...) pro- 
duces the elimination of the first-order character of the phase transitions; in other 
words, it eliminates the discontinuities in the thermodynamic expectation values of 
the conjugate quantities (magnetisation if the disorder is in the field, energy if the 
disorder is in the temperature, etc.). 

The problem was solved for the general case of spin variables a = {cTx} loc- 
ated on a D-dimensional lattice whose Hamiltonian is the sum of an ordered term 
(translation-invariant and non-random) and a fluctuating term with quenched ran- 
domness, represented in the following by a collection of independent random vari- 
ables {t]}. Some examples of this form are: 

1. Random field (RF) models 



where, in the ferromagnetic RFIM, cr G and J > 0. In the 0{N) model, ax 
are A^-component unit vectors with a rotation-invariant distribution. 




(3.1) 
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In RF models the spins are subjected to a fluctuating magnetic field composed 
of two terms: one uniform [h), and the other random, with the order of mag- 
nitude of e. We assume that the random fields rj^ are independently distributed 
with a probability measure u^drj) (with averages denoted by A{f) = f fu^drf)) 
that fulfil: 

A{r,) = Q , ^(r/2)>0, (3.2) 

and 

^(e"'') < oo , Vs < oo . (3.3) 



2. Random bond (RB) models 

For example, the Q-state Potts model, with a G {1, . . . , Q} and with Hamilto- 
nian with bond disorder 

Hl{,Cr) = ^{Jx-y + e:,-yVx,y)Sax,ay , (3.4) 

x,y 

or with site disorder 

x,y X y 

(3.5) 

3. Spin-glass models 

For example, the Ising model with Hamiltonian 

H{a) = 5Z ^^^y^^^v ~ + ^Vx)(^x ■ (3.6) 

\x-y\=l X 



In general, all the above models can be unified in a Hamiltonian of the form 

H{(t) = Ho{a) + J2 JL^f^- + ^aVa,x)9a{Tx(T) , (3.7) 
a X 

where the index a may parameterise pair-interaction terms of a given range or other 
multiple-spin terms, ga are bounded functions of the spin configuration, are 
translation operators (not to be confused with the temperature T = 1//3), and rja^x 
are a collection of random variables satisfying the conditions fl3.2p and fl3.3p . with 
an identical distribution within each a class. 

The free energy, F, is derived from the finite volume partition function Zy- By 
standard thermodynamic arguments, for almost every configuration of the disorder 
parameters {'r]a,x} the limit 

lim I \og[Zv{T, {h}, {e}, {r^})] = F(T, {h}, {e}) (3.8) 
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converges to a non-random function in the Thermodynamic Limit. In other words, 
the free energy self-averages. 

In addition, it is known that the free energy is convex in {h}, for fixed T and {e}, 
and therefore their directional derivatives exist; any discontinuity of those corres- 
ponds to a first-order phase transition. One can define the following order parameter: 



M,(T,{M,{e}) 



d 



d 



F{T,{h},{e}) 



In the case of the ferromagnetic RFIM 

1 



M(T,/i,e) = -[^((ao)+)-^((ao)-)] 



(3.9) 



(3.10) 



where with "-(-" and "— " we denote the extremal Gibbs states ("pure phases") 
constructed via choices of the boundary conditions (-1- or — ). 

The following is the main result for the general case, see Ref. |1UU] for its demon- 
stration: 

Theorem. In a D < 2 system with quenched disorder, described by a Hamiltonian of 
the general type of Eq. ( [5'. ?] ) with nearest neighbour interaction (this can be extended 
to longer range interactions) and with a continuous (non-atomic) probability measure 
u^drf) is 



M„(T,{/i},{e}) = 



V T > 0, {/i}, {e} and a for which > . (3.11) 



In the disordered Potts case, where the transition is due to a change in temper- 
ature, the free energy F is also convex in T for {/i} and {e} fixed and therefore its 
partial temperature derivative exists. We can define in this case the latent heat as 



(3.12) 



which is the order parameter for this first-order phase transition. In an analogous 
way it can be demonstrated that L = for D < 2 . 





( d 


d 






I 







3.2.2 Cardy-Jacobsen Theory 

The following results are based on a mapping between the Random Bond (RB) 
model (such as the dilute Potts model) and the Random Field (RF) model (such as 
the RFIM), see [lOlJ. Firstly we will summarise the main properties of the latter 
model. 

Random Field Ising Model (RFIM) 

It is defined by the Hamiltonian 

H = -j'^SiSj -\-'^hiSi -\- h'^Si , (3.13) 

i,j i i 
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hi being quenched random variables satisfying hi = and h"^ = A'^. The pdf of hi, 
p{hi), can be choserjfl Gaussian or bimodal (±/\). 

There are some important general theoretical results concerning this model: 

1. Dimensional reduction 

A D-dimensional system with a random field is equivalent to a system with 
D — 2 dimensions without the random field |lU8j . Therefore the lowest critical 
dimension is = 3 because for the pure Ising model it is Z)J."^ = 1. This 
result can be obtained by using supersymmetry arguments or through perturb- 
ation theories. Nevertheless, dimensional reduction seems to fail specifically 
for this model |109] . although it is a valid result for many other models. 

2. Imry-Ma argument 

Starting from T ^ in the RFIM and considering that the fundamental low 
temperature ferromagnetic state is Sj = +1, we can analyse what happens if 
we form a "droplet" with radius R with Sj = — 1 [IIOJ, see Fig. 13.11 



+ + + 




Figure 3.1: "Droplet" with different sign within an almost fully ordered Ising model. 
We study what is the effect of this perturbation depending on the dimension of the 
space. 

This "droplet" will present an interface with an energy cost 

AE = JR^-^ ~ R^'^ . (3.14) 
There is also an energy variation due to the random field within the "droplet" 

AEH = J2hi, (3.15) 

which, by the definition of the random field, will fulfil 

6=^J^ = ±iR''hlpy/^r^R^ , hlp = A\ (3.16) 

^Some controversy exists about the influence of this choice on the universality class of the model, 
see |107| and references therein. 
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We can always choose a point in the lattice where 5 < 0, so that the energy 
balance between AE and AEh produces the following results depending on 
the dimensionality of the space: 

• For D > 2 the fundamental low temperature state is stable. Consequently 
the low temperature ferromagnetic phase exists and a phase transition at 
finite temperature can be found. 

• For D < 2 the fundamental state is not stable and there will exist no 
phase transition. 

• For D = 2 we are in the marginal case. In Ref. [ITT] it was demonstrated 
that the rugosity of the interface (which is obviously not flat) destabil- 
ises the ferromagnetic state, see Fig. 13.21 This result is included in the 
Aizenman-Wehr theorem. 




Figure 3.2: Rough interface for a two-dimensional Ising model. 



The rugosity of the interface in D = 2 produces the energy [111] 

E = -C^R\ogR, (3.17) 
with C > 0, which makes the ferromagnetic state unstable. By defining (for D = 2) 

J{L) = JL- epilog L = JL{1 - CwlphgL) with: wip = ^ , (3.18) 

with L being the linear lattice size, we can obtain |109] 

dJ{L) 



dlogL 



J(L)(l-C<^) + 0«^), (3.19) 



which is the Renormalization Group (RG) equation for the coupling J. We can easily 
obtain the remaining RG equations by taking into account that h has dimensions of 
R^, h\p has dimensions of i?^, and J has dimensions of R^~^. Therefore it will be 
to leading order 
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( dh 



RF _ D 



dl 



^ h 



I ^ log(&) 



la 



(3.20) 



dJ_ 
, dl 



J[{D-l] 



(generalising to dimension D) 



In addition, by definition, wrf = h^p / J and therefore 

dWRF 



dl 



-Wrf + Cwjip , witli: e = D — 2 . 



(3.21) 



If D > 2 then e > and the RG flow has a non-trivial random fixed point (with 
w ~ 62 ^ 0) in agreement with [llUj . Using these results, the phase diagram for 
the RFIM can be obtained, see Fig. 13.31 




Figure 3.3: RG flows for the RFIM for D > 2 . 



Cardy-Jacobsen mapping 

In the case of a pure system undergoing a first-order phase transition there will be 
coexistence of a (generally unique) disordered phase and the (generally non-unique) 
ordered phases. The internal energies Ui and U2 of these two phases differ by the 
latent heat. Consider, see Ref. [lOlj . a large (say horizontal) interface between the 
disordered phase and one of the ordered ones, with surface tension E. If 3> 1, 



59 



Chapter 3. Quenched Disorder Effect on a First-Order Phase Transition 



there will be a really small number of isolated bubbles of the opposite phase above 
or below the main interface. The free energy of these bubbles is proportional to 
their areas multiplied by the surface tension. 

On the other hand, let us consider an Ising model and build an interface between 
the two possible ordered phases (all spins taking the values ±1). At a very low 
temperature, there will basically exist no "bubbles" above or below the interface and 
the surface tension will be ~ 2 J, where J is the reduced coupling of the Hamiltonian. 
In the limit ~ 2 J ^ 1, these two interface models will be identical. 

We will analyse the effect of disorder in these two models. In the first one, we 
introduce random bonds (disorder coupled with the energy) while in the second we 
introduce a random field (disorder coupled with the magnetisation). The changes 
in the energy due to the introduction of the disorder are: 

1. Random Field: 

5^ Mr) Mr), (3.22) 

where the sums are defined over all the points above (>) or below (<) the 
interface. 

2. Random bonds 

Ui Sx{r) + t/2 5^ 5x{r) = Ij^ ^^i^) - ^ 5x{r) j + const. , (3.23) 

with 6x{r) being the local impurities density, L = Ui — U2 is the latent heat, 
and the final constant is md^pendent of the interface locafcfl The latter 
equation has the same form as that corresponding to the random field. 

Therefore the thermal variables of the random bond system are related to the 
magnetic variables of the RFIM by the following mapping: 

Random Bond Random Field 

E/kT^ J/kT 

{L/kT,)x ^ hRp/kT (3.24) 
(T - TJL HM 



The energy can be split in the following way, A being the difference between the two sides of 
Eq. (1^:^ : 

^=IU2Y^ 5x{r) + \u, 5x{r) + ^ 5x{r) + ^ 5x{r) = 



dx{r) + 6x{r) 



Y^5x{r)+Y.^<'') 



^ 5xi + i[/2 ^ 5xi 



which is independent of the interface location. 
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The last relationship is between the "field" (T — T^) and a uniform external field 
H, which helps to distinguish between the two phases (in the same way as [T — T^). 
One of the possible problems of this mapping is its use of the local energy density as 
a kind of order parameter. However it can be made completely explicit, for example 
for the Q-state Potts model, through the mapping to the random cluster model 
where ~ L ~ logQ, Q ^ oo, see Ref. [lUlj . 



Explicit relationship with the Potts model 

We can now derive |101j the specific relationship of the previous section with the 
Q-state Potts model with quenched disorder, with Hamiltonian 

H = - ^ Ki,5s,s, , (3.25) 

<i,j> 

where the sum extends only over nearest neighbours. The ferromagnetic couplings 
Kij are quenched random variables, taking the values Ki and K2, each with prob- 
ability |; in other words, their pdf is 

p{K,,) = ^6iK,, - K^) - ^6{K,, - K2) . (3.26) 

When (e^^ — l){e^'^ — 1) = 1 this model is, on average, self-dual, and, if the 
transition is unique, is at its critical point |112] . It is useful to parameterise the 
model through 

u^J = (e""- - 1) = Q^+"-- , (3.27) 

with Wij = ±w. w > measures the strength of the randomness, with w = being 
the case without disorder. We can solve for 

ir,, =log(l + Q^+"'-), (3.28) 
and consider the limit Q — )■ 00 to approximate 

= log(l + qH"--) ~ log(Q5+->.) = Q + nj,,^ logQ . (3.29) 
By substituting in the Hamiltonian we obtain 

^ = l0g<5 ^s,s, - logQ ^ W^j6s,s, , (3.30) 

<i,j> <i,j> 

where the first term corresponds to an ordered model while the second term corres- 
ponds to a disordered one. Therefore the term added to the pure model is 

Padded ~ ^ Wij log Q 6s,sj ■ (3.31) 
<i,j> 

We will work in the following in two dimensions; i.e. the label i identifying the 
lattice sites means i = {x,y), with {x,y) being Cartesian coordinates. Using the 
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same arguments as in the previous section, if an interface divides the space into two 
parts, denoted by ">" and "<", then 

Padded ~ ^wlogQ + ^wlogQ, (3.32) 

> < 

noting that within each of these two homogeneous regions Sg-g. = 1. Each of the 
two terms of Eq. fl3.32p is 

E»"'SQ = 5(E + E)+5(E-E). (3-33) 

but, as in the previous section, the term {J2> + J2<) (^oes not depend on the inter- 
face position and as a consequence 

^(E-E]^l°g^ (3.34) 

is an analogue of the random field term of Eq. (13.221) . Therefore we obtain the 
relationship 

hRF < — > -w log Q . (3.35) 
In addition, following Ref. [S2], when Q ^ oo, the surface tension is: 

r~llogQ. (3.36) 

But as was seen in Sec. 13.2.21 for the RFIM ~ 2 J, and therefore we have the 
relation 

J^^logQ. (3.37) 

Finally, while a uniform field in the RF model distinguishes between the two 
phases, in the RB model this is the task of the reduced temperature t = 
provided that t is coupled to the energy density we can make the identification 

hi — ^^tlogQ. (3.38) 
To summarise, the mapping between the Potts model and the RFIM is: 



RFIM Q-state Potts model 

J ^ |logQ 

h ^ i^logQ (3.39) 

flRF log Q 



We can use this mapping to derive the RG equations for the Potts model starting 
from those for the RFIM [109j , see Eq. f l3.20p . We will denote for the sake of clarity 
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the w used in Eq. fl3.35p as w^ib = w, to distinguish between this w^b and the w^p 
defined in Eq. fl3.18p . Therefore from Eqs. fl3.35p and (13.370 one obtains 



But wjip = hjip/J, hence 



wrb = -^wrf , (3.41) 
and therefore the RG equation for wrf is also vahd for wrb, i-e., 

= -^WRB + AwIb , with: e = D-2 , 1 = log(6) . (3.42) 

From the last RG equation of Eq. (I3.20p and from Eq. (I3.37p . one easily derives the 
RG relation 



-(logg)-1p-l)-A<B]- (3.43) 
Finally, using again Eq. (I3.20p and Eq. (I3.38p . one obtains: 

I = t(l + AwIb) ■ (3.44) 

Summarising, the set of RG equations for the Q-state Potts model with quenched 
disorder is: 



r dw 



RB 



dl 



-^wrb + Aw 



\q , with t = D — 2 



I ^ \og{h) 



dt 



dl 
dl 



t(l + Awl 



RBl ' 



(3.45) 



-{\ogQ)-^[e + l-AwlB\ 



We can now analyse the stability of the fixed points of the RG transformations 
to trace the phase diagram of the model. To simplify we make the change 



logg 



Then the fixed points of the RG transformation are the Gaussian one: 

Q = t = Wrb = , 

and the tricritical one: 



Q = t = 



Wrb 




(3.46) 



(3.47) 



(3.48) 



63 



Chapter 3. Quenched Disorder Effect on a First-Order Phase Transition 



The Jacobian matrix for the transformation is 

/-f + 3Aw;2„ 



J 



\ 2AQwrb 
which evaluated at the Gaussian point is 



2AtwRB 1 + Aw%B 

-(e + l-Awl^)/ 



J\ 



Gauss 



/-f 

1 








\0 -(e+1)/ 

resulting in that at the Gaussian fixed point the "fields" wrb and Q are irrelevant 
while t is relevant. At the tricritical (TC) point, the Jacobian will be 



TC 







1 + f 



\0 



\ 



;i + f)/ 



and therefore wrb and t are relevant "fields" while Q is irrelevant. For D > 2 the 
eigenvalues are 



Vw 

yt - 
vq = 



e > 



1 + - 



Relevant 
Relevant 

Irrelevant 



The {Q,w) 
system, for Q 



plane of the phase diagram is depicted in Fig. (13. 4p . In the pure 
> Q2 there is a phase transition with non-vanishing latent heat 
controlled by a fixed point at infinite Q. For D > 2 it will continue into the 
shaded region bounded by a line of tricritical points whose exponents are related 
to those of the RFIM [101] . It also may be shown that the latent heat vanishes as 



w 



as the line RQ2 is approached from below. Let us note also that Q 



corresponds to the percolation model, where the disorder is irrelevant (a ~ —0.64), 
while Q = 2 corresponds to an Ising model, where the phase transition is always of 
second order but the disorder is relevant {a ~ 0.1118). Therefore there exists a Qi, 
with 1 < Qi < 2, at which the sign of a changes. In addition, for Q > Q2 = 2 + e 
the transition becomes first order. 

Above the line RQ2, as w grows the RG equations lose validity. In addition the 
surface tension goes to zero and the mapping between the two models disappears. 
Nevertheless, for Q ^ 00 the mapping remains exact and the flow goes to infinite w. 
But this can not happen for finite Q because this is the percolation limit K1/K2 = 0, 
at which the disorder is relevant |113] . Therefore there should exist |101] another 
line of stable fixed points emerging from Pi, which control the universal continuous 
transition for large, but finite, values of w and Q. 
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Figure 3.4: Phase diagram of the dilute Potts model with Q states for D > 2, 
obtained in Ref. [lUlj . Q grows towards the left; w is a measure of the disorder 
strength, w = being the pure case. The latent heat do not vanishes within the 
shaded region while outside the transition is continuous. Pi and P2 are the percol- 
ation thresholds. 



3.3 The Model 

In the three-dimensional site-diluted Q-state Potts model [H2] the spins, cTj, take the 
values (Tj = 1, . . . , Q and are defined at the nodes of a cubic lattice with probability p. 
We consider only nearest neighbour interaction and periodic boundary conditions. 
Therefore the Hamiltonian takes the form: 

•H^P'" = - J]e,e,<5^,.,, (3.49) 

with ej being quenched occupation variables (e, = or 1 with probability 1 — p and 
p respectively)^!, and denoting nearest neighbours. Each one of the specific 

disorder realizations ({e} spatial distribution) is called a sample. The pure system 
is recovered for p = 1, and is known to undergo a first-order phase transition for 
Q > 3 [inilin] generally regarded as very strong. 

A valid order parameter for the model is the magnetisation density (a Q-dimensional 
vector) defined as 



^To reduce statistical fluctuations, we kept only the spins in the percolating cluster |114| that 
control the critical behaviour. However, in the most interesting region (p > 0.9) this correction is 
quite small. 



QS 



(3.50) 
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with V = being the volume and L the hnear size of the system. We can define 
the magnetic susceptibihty as 

X = V\M\\ (3.51) 

A well-behaved definition for the correlation length in a finite system is obtained 
from the correlation function as [7] 

^ ^ (3.52) 



^4sin^(7r/L) 
where 

V- 



F ^ -(|F(27r/L, 0, 0)|2 + |F(0, 27r/L, 0)|2 + |F(0, 0, 27r/L)|2) , (3.53) 

and where we denote the thermal averages with brackets while the sample average 
is overlined. In addition 

F{k) = ^J2^""''^rcrr. (3.54) 

r 

In this work we use the microcanonical simulation method defined in Ref. |13] . 
see also Sec. I2.2.2[ so that, by using the Maxwell construction, see Appendix [El we 
can directly obtain some quantities characteristic of the phase transition. Firstly 
the critical temperature is fixed by the definition of the Maxwell construction: the 
e-integral of /9{e}(e) — 1/Tc from to Cq must vanish, where e is the energy density 
and /3{e}(e) was defined in Eq. (12. 58 p . This fact also implies that 

Sd - So = (ed - eo)/Tc , (3.55) 

s being the entropy density. The latent heat is defined directly as 

Ae = ed — Co . (3.56) 

Finally the surface-tension, E, is L^/2 times the integral of the positive part of 
/3{,}(e) - l/Tc, see Ref. [I3]. 



3.4 Numerical Results 

We have studied numerically two cases of the three-dimensional site-diluted Potts 
model: the four-state {Q = A) and the eight-state {Q = 8) cases. Both cases undergo 
a well-known |82j strong first-order phase transition in the pure case; the strength 
of the first-order character of the phase transition will grow with Q. In both cases a 
softening of the discontinuities of the first derivatives of the free energy is expected 
to appear with increasing dilution. The critical concentration, pc, at which the 
character of the phase transition switches to second order will depend on Q, being 
smaller for increasing Q. 
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3.4.1 Methods 
Simulation method 

To update the Potts spins of our systems, we used a microcanonical version [13] 
of the Swendsen-Wang (SW) [llj cluster method. For disordered systems, SW up- 
dates loosely connected regions properly |ST] and does not require tedious parameter 
tuning. The microcanonical cluster method, which is not rejection- free, depends on 
a tunable parameter, k. In order to maximise the acceptance of the SW attempt 
(SWA), K should be chosen as close as possible to /3{e}(e)- After every e change, we 
performed cycles consisting of 10^ Metropolis steps, a k refresh, then 10^ SWAs, and 
a new k, refresh. The cycling was stopped, and k fixed, when the SWA acceptance 
exceeded 60%. We then performed a number of SWAs depending on the lattice size, 
taking measurements every 2 SWAs. 

For Q = 4 we performed thermalization checks that included comparisons of hot 
and cold starts or even mixed configurations (bands and strips [13] )• We checked that 
the Maxwell construction obtained for the pure case of the largest system (L = 128) 
does not depend on the initial configuration, after discarding a part of the initial 
Monte Carlo history. 

In the Q = 8 case, reaching the thermodynamic equilibrium is a far more com- 
plicated task, especially on the first-order side of the phase diagram. For a first-order 
phase transition, it is known that metastable states do exist. These states can have 
a very long life even when they are not the true equilibrated states. This is most 
dramatic for large systems in which the simulation times are intrinsically longer. 
Therefore the thermalization issue in this case deserves a special treatment that will 
be described in Sec. 13.4.31 

Sample averaging methods 

For a disordered system, one has to analyse a set of functions /3{e}(e) corresponding 
to a large enough number of samples. There are two natural possibilities. On the 
one hand, one can use the Maxwell construction for each sample extracting T^, e^, 
Co and S, and then considering their sample average, median, or even their pdf, see 
Fig. 13.61 This is the most usual approach. 

On the other hand, one can compute the sample average of the (inverse) tem- 
perature defined at each simulation energy, e, /3(e) = /3{e}(e), and then perform the 
Maxwell construction on it (i.e., take the sample average of s(e), rather than the 
sample average of the free energy at fixed T). 

We found empirically that the two sample averagings are equivalent in the first- 
order piece of the critical line. This is hardly surprising, because the internal energy 
as a function of T is a self-averaging quantity, for all temperatures but the critical 
one. Therefore, also Cd, Co, and Tc are self- averaging properties in the first-order 
part of the critical line. 

While the first method offers more information, it is computationally more de- 
manding (it requires high accuracy for each sample). The method featuring /3(e) 
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can be used as well in the second-order part of the critical line, but its merit in that 
region is yet to be investigated. 



3.4.2 D = 3, Q = 4 Site-Diluted Potts Model 

For Q = 4 we investigated the phase transition for several p values in the range 
0.75 < p < 1. As a rule, we found that at fixed p the latent heat is a monotonically 
decreasing function of L, see Fig. 13.71 For each p value, we simulated L = 16, 32, 
64, and 128 (for a given p, we did not consider larger lattices once the latent heat 
vanished). For all pairs {p, L) we simulated 128 samples. Also, some intermediate p 
values were added for the FSS study, see Fig. 13.81 and we raised to 512 the number 
of samples for (L = 16 and 32, p = 0.86 and 0.875). 



General behaviour 

Following the procedure of Sec. 13.4.11 we performed sample averages of the Maxwell 
constructions for each L and p. In Fig. 13.51 one can see the general behaviour of 
the Maxwell constructions as the spin concentration, p, varies: while for large p 
{p ^ 0.95) we can form the Maxwell construction for every system size, it softens 
with decreasing p up to a point, Pc, at which both latent heat and surface tension 
vanish. This Pc depends on the system size. 
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Figure 3.5: Sample-averaged e-derivative of the entropy, /3(e), for several lattice 
sizes, L, and spin concentrations, p. Metastability requires a non- decreasing /3(e). 
The horizontal line marks the critical (inverse) temperature 1/Tc, obtained through 
Maxwell construction. At fixed L the surface tension increases for increasing p. Note 
that, for a fixed dilution, a seemingly first-order transition [L = 64, bottom-right), 
may actually be of second order if studied on larger lattices {L = 128, bottom-left). 
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Figure 3.6: Histograms for the sample- dependent latent heat Aj^je = Cd — Cq (left) 
and surface tension, Z'(right). In the top panels we show results in the largest 
lattice, where two very close spin concentrations behave very differently. The three 
types of horizontal lines drawn (indicating central value and statistical error) cor- 
respond, from top to bottom, to the median, the mean, and the value obtained from 
/3(e). In the lower panels we show the histograms for p = 0.98 and different L's 
(note the difference in the horizontal scales with the upper part). As can be seen, 
the latent heat is self- averaging while the surface tension is not. 



As was said before, for each sample we can define the different thermal- averaged 
quantities, and then determine their mean, median, or pdf. We can also compute 
the sample- average /3(e) = /3{e}(e), and then perform the Maxwell construction on 
it. We compared the two approaches for this model both for the latent heat and for 
the surface tension, see Fig. 13.61 We found that the two approaches are equivalent, 
although the second one requires less statistical accuracy for each sample and is 
therefore less numerically demanding. Also from Fig. 13.61 (top row) we can see 
that as the dilution is slightly decreased (the tricritical point is reached), a great 
number of samples present vanishing latent heat and surface tension; the transition 
has become continuous. Finally we found that within the first-order part of the 
phase diagram {p = 0.98) the width of the histograms of the latent heat decreases 
as the lattice size increases; this is the definition of a self-averaging quantity. On 
the contrary, we can not see this behaviour for the surface tension and therefore we 
can state that it is not self-averaging. 
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Figure 3.7: Top: Latent heat as obtained from /3{e) as a function of spin concen- 
tration for several lattice sizes (lines are linear interpolations). Data for p = 1 and 
L = 128 were taken from Ref. [13]. To illustrate the sample dispersion, we also show 
the scatter-plot of {N/L^, A{e}e) for the 128 samples at L = 16 j9 = 0.85 and L = 64 
j9 = 0.92. Bottom: As the top panel, but for the surface tension. 



Latent heat and surface tension 

Our results for the behaviour of the latent heat and the surface tension obtained 
from /S(e) as the dilution changes are shown in Fig. 13.71 The apparent location 
of the tricritical point (i.e., the p where both Ae and E vanish) shifts to higher 
p for increasing L rather fast. For lattice sizes comparable with those of previous 
work |16], L = 16, we obtain a sizeable value p^^^^ ~ 0.75, but the estimate of pc 
increases very rapidly with L. An extrapolation to L — )■ oo is called for. 

The pdf 's for Ae and E, Fig. 13. 6^ display an interesting L evolution. When 
the /3(e) changes behaviour from non-monotonic [L = 64, Fig. 13. 5^ bottom-right) 
to monotonic (L = 128, Fig. 13. 5[ bottom- left), the two pdf's becomes enormously 
wid^, see the top panels in Fig 13.61 This arises because for many L = 128 samples 
the curve /9{e}(e) is becoming fiat, or even monotonically decreasing (i.e., Ae = S = 
0), while no such behaviour was seen for L = 64. Only for p = 0.98 does the width of 
the pdf's for Ae scale as as expected for a self- averaging quantity, see Fig. 13.61 

- bottom-left. The surface-tension is not self-averaging, see Fig. 13.61 - bottom- right. 

Finite Size Scaling study 

From Figs. 13. 5| 13. 6[ and 13. 71 one cannot rule out that Pc ^ a disordered first-order 
transition would not exist. Fortunately we can solve this dilemma by consider- 

^The estimates for Ae and U are consistent with the median of their (non-Gaussian) pdf's. 
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Figure 3.8: Correlation length in units of the lattice size in the Q = 4 case, at phase 
coexistence for the paramagnetic phases, function of concentration, p, for 

several system sizes, L (lines are cubic spline interpolations for data at fixed L). 



ing the correlation length, obtained from the sample-averaged correlation function, 

Eq. (issa). 

We take the correlation length in units of the lattice size at Cd (see Fig. 13. Sp . and 
Co (see Fig. 13. 9p . as obtained from /3(e) (a jack-knife method [7j takes care of the 
statistical correlations). For all p < pc, one expects that both ^{cdj/L and ^{eo)/L 
tend to non- vanishing and different limits for large I^. For p > pc, $,{ed)/L is of 
order 1/L, while ^{eo)/L ~ L^/^. For a fixed L, with increasing p, the behaviour 
goes from a second-order type to first-order (see Fig l3.5l) . Hence, a FSS approach [7j 
is needed. 

Consider the curves of ^{e^)/L versus p for different L, see Fig. 13.81 There is a 
unique concentration, p^'^^^ where the correlation lengths in units of the lattice size 
coincide for pairs of lattices of sizes L and 2L. One has |fl 

^ p, + AdL-^ (3.57) 

An exactly analogous result holds for ^ (cq) / see Fig. 13. 9[ Since and are quite 
different, see Fig. I3.10[ a combined fit of all the data yields an accurate estimate of 

^We have checked numerically that this is indeed the case for the D ^ 2, Q — A, pure Potts 
model (a prototypical example of a second-order phase transition with a double-peaked canonical 
pdf for e at Tc), see Sec. [2X3l 

®The tricritical point has no basin of attraction for the RG flow in the (T,p) plane. Although 
two relevant scaling fields are to be expected, the Maxwell construction allows us to eliminate one 
of them and hence we use the formulae for a standard critical point. 
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Figure 3.9: Correlation length in units of the lattice size in the Q = 4 case, at phase 
coexistence for the ferromagnetic phases, Cq, as a function of concentration, p, for 
several system sizes, L. 




Figure 3.10: Spin concentration where the values of ^/L (data from Figs. 13.81 and 
13. 9p coincide for lattices L and 2L versus l/L^, see Eqs. f l3.57p and (13. 58 p . Lines are 
a combined fit for x, Pc, A^, and A^. 
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the location of the tricritical point: 

Pe = 0.954(3), X = 1.23(9), = — -, C.L. = 24% . (3.58) 

Of course, due to higher-order scahng corrections, Eq. fl3.57p should be used only for 
lattices larger than some L™'" [70]. The fit 'was acceptable taking = 12 and 
L™° = 16 (for the sake of clarity we do not display data for L = 12 in the figures). 
Therefore we can conclude that p = 0.98 is definitively in the first-order part of the 
critical line. 

We now look at at p^''^^^ see Figs. l3.8l and [Xin Consider ^ (cd) /L as a function 
of {L,p), see Fig. 13.111 Its salient features are: 

1. For fixed L, ^{cdj/L is a decreasing function of p (while ^{eo)/L is increasing). 

2. For fixed p, ^{cai/L has a minimum (while ^{eo)/L has a maximum) at a 
crossover length scale, L^rip), that separates the first-order type of behaviour 
from the second-order type, see Figs. 13.81 and 13.91 

3. At the crossing point p^^"^^ we have L < Lcr{p^'^^) < 2L. 

4. At least within the range of our simulations, L^rip) is an increasing function 
of p. 

A standard scaling argument, combined with (1-4), yields that C,{ed)/L at p^'"^^ is 
of order 1/L„(p^'2^) {^{eo)/L ~ LcJ"^). If L„(p) diverges at Pe, ^{ed)/L at p^'^^ 
should tend to zero for large L, which is indeed consistent with our data. 

3.4.3 D = 3, Q = 8 Site-Diluted Potts Model 

We also present in this chapter some of the preliminary results of our study of the 
eight-state {Q = 8) site-diluted Potts model using basically the same methodological 
approach as in the Q = 4 case. First, however, it has to be stressed that there are 
two important differences in this case: 

• We used chiefly another kind of computing platform. While the Q = 4 case 
was entirely simulated on typical cluster facilities, i.e., the BSC and BIFI, 
the Q = 8 case case was simulated on IBERCIVIS, a distributed computing 
platform based on BOINC, see Appendix O This change in platform involved 
both advantages and disadvantages. By using IBERCIVIS, we were able to 
outperform broadly all previous statistical accuracies both in the number of 
samples (we were able to simulate up to 2000 samples of a system with 64^ 
spins) and in the number of dilution levels (around ten for each system size). 
In addition, we did more than 3 x 10® Swendsen-Wang steps at each energy 
of a system with 64'^ spins. In particular, with IBERCIVIS we obtained more 
than 300 years of computation time in less than a year of wall clock time. This 
would have been hard to achieve using a traditional cluster facility. 
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Figure 3.11: Correlation length at in units of the lattice size, for fixed dilutions 
as a function of the inverse lattice size. For fixed p < pc, C,{ed)/L has a minimum 
at a crossover length scale, Lcr{p), that separates the first-order type of behaviour 
from the second-order type. 




On the other hand, the use of IBERCIVIS, a novel infrastructure, led to nu- 
merous unusual problems in the adaptation and stabilisation of the original 
code to the new computing paradigm, see again Appendix [Gl The huge output 
of the computations has to be carefully analysed, and a major effort must be 
made to identify all the possible error sources. For example, the connection 
of the individual parts of each BOINC job is an extremely delicate issue, and 
the development of a secure mechanism for the detection of corrupted outputs 
is fundamental. 

• It is known that the pure Potts model undergoes a first-order phase transition 
for Q > 2 in three dimensions [82] , with the strength of the first-order character 
being larger as Q grows. Therefore the Q = 8 case will show more evidently 
the features of this kind of transition (i.e., latent heat, metastabilities, phase 
coexistence, etc.). This fact has pros and cons. The main benefit is that if 
we want to see a first-order phase transition in the presence of disorder, the 
first-order region in the phase diagram is expected to be larger for Q = 8 than 
for Q = 4, in other words Pc will be smaller. This will mean stronger evidence 
for the main result of this chapter - first-order phase transitions do exist in 
the presence of disorder in three-dimensional systems. 

On the contrary, the stronger first-order character of the transition produces 
much more palpable metastability effects. This is a huge problem in Monte 
Carlo simulations, because exponential autocorrelation times will grow sub- 
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Figure 3.12: Comparison of the different "energy walks" for the system with L = 48 
and p = 0.95. The upper part correspond to the Maxwell construction of a random 
sample. Note the clear difference between the walk starting from hot (red solid line) 
and cold (blue dashed line). The lower part shows the softening of the difference 
as the sample average (with 40 samples) is performed. Horizontal lines mark the 
corresponding estimates of (3c- 



stantially, see Appendix O and Ref. |115j . making thermalization really hard 
to achieve for large systems at given values of their internal energy. This fact 
restricted us to thermalizing systems with "only" 64'^ spins on the first-order 
side of the phase diagram, compared with the Q = 4 case for which we were 
able to thermalize systems with up to 128^ spins. Anyway we plan to study 
the first-order side of the transition by making estimates of the errors due to 
not having reached the asymptotic states of the system. 

For Q = 8 we simulated the model for p values in the range 0.65 < p < I. For 
each p value, we simulated L = 12, 16, 24, 32, 48, 64, and 96 (for a given p, we did not 
consider larger lattices once the latent heat vanished). Finally, we disregarded our 
simulations for L = 96 because of the impossibility of thermalization in a reasonable 
time. For all pairs (p, L) we simulated at least 500 samples. 

To check the thermalization of the systems we compared simulations of the same 
samples (distribution of the vacancies), performing annealings starting from both 
random configurations at high temperatures and cold (all the spins in the same state) 
configurations at low temperatures. As we performed our "energy walk" we found 
that for energies corresponding to pure states (with no "islands" of the other phase) 
both annealings will agree fully. However, between Co and Cd there will exist some 
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energies where the two anneahngs will produce different estimates of the observables 
(especially for the largest lattices). These energies are precisely those at which the 
system switches between the different configurations of the "islands", for example 
from a "droplet" to a "strip" . In the case that the thermal averages of the different 
observables from the two annealings were similar, we would be fairly confident of 
the equilibration of the system. If they were not, we could at least estimate the 
error due to the lack of thermalization from their difference. 

In Fig. 13.121 we plot the comparison of the annealings of the system with 48^ 
spins and p = 0.95. The system is clearly undergoing a first-order phase transition. 
The simulation of each sample used in Fig. 13.121 took around three days of a last 
generation Pentium 17 3.0 GHz core and are clearly not thermalized! With this in 
mind, it is clear that on the first-order side of the phase diagram the simulations 
must be really long in time to reach equilibrated states. As was said before, this 
fact will critically restrict us in simulating large systems. 

Due to the lack of mixed phases, thermalization is quite easy to achieve on the 
second-order side of the phase diagram. No metastability will exist and the cluster 
update method will work very well. 

Therefore the approach to the problem must be very different depending on the 
dilution of the system. If the dilution is weak, the systems will undergo first-order 
phase transitions and we will not be able to simulate large systems. Nevertheless we 
can estimate the latent heat and the surface tension to obtain the exponents of their 
scaling, always taking into account the possibility of unequilibrated systems. On the 
contrary, if the dilution is important, since the systems will undergo second-order 
phase transitions, we will be able to equilibrate large systems (with 64^ spins) and 
to obtain accurate results for the tricritical point location. In this work we present 
only the latter study, i.e., the study of the exact location of the tricritical point. 
The study of the first-order side will be left for further research. 

General behaviour 

First we outline the behaviour of the model, which is very similar to that of the 
Q = 4 case although the first-order character is stronger. Firstly, by taking the 
sample average of the Maxwell constructions, see Sec. 13.4.11 we can obtain for each 
system size the behaviour of the model as the dilution changes, see Fig. l3.13l for the 
system with L = 24. Note that in this case we obtain clean Maxwell constructions up 
to p = 0.800. The system is undergoing (on average) a first-order phase transition 
even with 20% of vacancies! It is also remarkable that in the pure case, p = 1, 
the Maxwell construction is smooth, without flat parts or strong steps between 
consecutive energies. 

By representing the same plot for a larger system, see Fig. 13.131 for the L = 48 
case, we find that Maxwell constructions can not be formed for dilutions less than 
p = 0.850. As was said in Sec. 13. 4. 2^ the apparent tricritical point depends on the 
system size (a FSS study is again called for). We can also see that the Maxwell 
construction in the pure case presents fiat parts and clear steps in the temperature. 
They are due to the existence of clear mixed regimes with droplet or strip-like 
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Figure 3.13: Maxwell constructions for the system with L = 24 as a function of the 
spin concentration, p. In the dilute cases, each point represents the average of the 
temperature for 500 samples. Note that for this size we can form the construction 
even for p = 0.800. 




Figure 3.14: Maxwell constructions for the system with L = 48 as a function of the 
spin concentration, p. In the dilute cases, each point represents the average of the 
temperature for 500 samples. Note that we can form the construction only up to 
p = 0.850, the transition has become of first order between p = 0.850 and p = 0.800. 
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Figure 3.15: Top: Latent heat as obtained from the averaged Maxwell construction 
as a function of p for each simulated lattice size (lines are linear interpolations) for 
the (5 = 8 case. Bottom: As in the top panel, but for the surface tension. 



configurations in which the internal energy basically does not change over a range of 
energy densities. This is also manifest in the case of the mildly dilute samples prior 
to the sample- averaging process. The location of the fiat parts and the steps is by 
far the part of the "energy walk" that is most sensitive to the metastability effects; 
it is really difficult for any Monte Carlo spin update method to perform properly 
with configurations of this kind. 

Again we can plot the entire behaviour of the latent heat and the surface tension 
as a function of the system size and the dilution, see Fig. 13.151 This figure must be 
compared with that corresponding to the (5 = 4 case. Fig. 13.71 While in the (5 = 4 
case the tricritical dilution is clearly above p = 0.95, this is not so for (5 = 8. The 
data for both the latent heat and the surface tension show the clear trend of the 
tricritical dilution towards larger values as the number of Potts states grows. The 
points in Fig. 13.151 corresponding to large values of the dilution and the lattice size 
[j) > 0.925, L > 32) are possibly not fully equilibrated, so special treatment of the 
data is needed to obtain accurate information of the scaling in this part of the phase 
diagram. This will be done in future work. 
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Figure 3.16: Correlation length for the Q = 8 case in units of the lattice size, at 
phase coexistence for the paramagnetic phases, function of concentration 

for several system sizes, L. The huge error bar for L = 48, p = 0.85 is due to 
the enormous sample-to-sample dispersion of this observable close to the tricritical 
point. 

Behaviour of the model on the second-order side 

From Fig. I3.15[ while one can not obtain an accurate estimate of pc ^ 1, one can 
again consider the correlation length obtained from the sample- averaged correlation 
function, Eq. fl3.52p . We use then the same approach as we used in Sec. 13.4.21 
computing the crossings of the correlation length in units of the lattice size at Cd, 
see Fig. I3.16[ and Co, see Fig. I3.17[ as obtained from /3(e). 

As was done in Sec. I3.4.2[ we define p^^^^ as the crossing points of the correlation 
length ^(cd) (in lattice size units) for pairs of lattices with L and 2L. We can fit 
these points again to the form of Eq. (I3.57P to obtain the value of the tricritical 
dilution pc. Fitting our data set in Cd, see Fig. 13. 16^ we obtain: 

Peed = 0.915(1), x = 1.14(28), ^^ = ^> C.L. = 47%, (3.59) 

which is a perfectly valid fit producing a value for p^ clearly less than unity. Using 
the same approach for the crossings of ^(co), see Fig. 13.171 we obtain a valid fit with 
parameters 

y2 1.1 

Pceo = 0.910(2), x = 0.95(49), = C.L. = 58% . (3.60) 

Finally we can fit both data series to the form f l3.57p sharing the same coefficients 
Pc and X. To get an acceptable value for the of the fit we had to disregard the 
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Figure 3.17: Correlation length for the Q = 8 case in units of the lattice size, at 
phase coexistence for the ferromagnetic phases, Co, as a function of concentration 
for several system sizes, L. 



data with L < 16 for ^(cd) and L < 24 for .^(eo) obtaining the fitting parameters: 

1.50 

Pcjoined = 0.922(1), Xjoined = 0.93(47), ^ = C.L. = 47% . (3.61) 

A plot of all the above fits is shown in Fig. 13.181 Therefore we can firmly conclude 
that p = 0.925 is in the first-order part of the critical line for the three-dimensional 
Potts model with Q = 8. This is another result that reinforces our main conclusion 
of the previous section, i.e., first-order phase transitions do exist in D = 3. 



3.5 Conclusions 

In this chapter we have performed a detailed study of the effects of quenched dis- 
order on a three-dimensional system undergoing a first-order transition in the pure 
case. We studied the site-diluted version of both the Q = 4 and the Q = 8 Potts 
model, a model undergoing a prototypically strong first-order transition, with the 
strength being proportional to the value of Q. A small degree of dilution smooths 
the transition up to the point of becoming second order at a tricritical point, p^. 
We observed strong finite-size effects in both the location of the tricritical point 
and the behaviour of the most relevant quantities (latent heat, surface tension, cor- 
relation length, etc.). A dehcate FSS analysis allowed us to firmly conclude that 
Pc < 1, with pc = 0.954(3) and pc = 0.922(1) in the Q = 4 and Q = 8 cases respect- 
ively. We are then able to claim that (quenched) disordered first-order transitions 
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Figure 3.18: Spin concentration where values of ^/L (data from Figs. 13. 1^ and IXT7|) 
coincide for lattices L and 2L versus 1 /L^'j™°<='^ , see Eq. fl3.6ip . Solid lines are the 
combined fit, while the dashed and dotted curves correspond to the individual fits 
of the crossing of ^(cd) and ^(co) respectively, Eqs. (13.591) and (13.601) . 



do exist in three dimensions, although quenched disorder is unreasonably effective in 
smoothing the transition (we speculate that the percolation mechanism for colossal 
magnetoresistance proposed in [21] could be fairly common in D = 3). 

We also observed that, for a given p < Pc, a crossover length scale Lcr{p) exists 
such that for L < Lcrip) the behaviour is of first-order type. The asymptotic second- 
order behaviour appears only for L > L^rip). 

In the Q = 4: case, we also verified that the latent heat is a self- averaging quantity 
for random first-order phase transitions while the surface tension is not. We will try 
to verify this point for the Q = 8 case in future work. 

All these results were made possible first by a new definition of the quenched 
average that avoids long-tailed pdf's [?6J, and second by the use of a recently intro- 
duced microcanonical Monte Carlo method that features the entropy density rather 
than the free energy [T3] . 

As further research, we will obtain novel information on the scaling of some 
quantities on the first-order part of the critical line in the (5 = 8 case. To perform 
this analysis, we will have to deal with systems that are not fully equilibrated. The 
characterisation of the effects due to the metastable states will be done by comparing 
pairs of simulations performing annealings from hot and cold states. 
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The Site-Diluted Heisenberg 
Model in Three Dimensions 

4.1 Introduction 

The three-dimensional Heisenberg model is the most general representation of the 
interaction of the spins within an isotropic magnetic material, where isotropic means 
that the magnetisation does not have any preferential direction to point to. Besides, 
other popular models such as the Ising or XY models describe materials with a 
plane or axis of easy magnetisation, as is the case for instance of hexagonal lattices 
where the magnetisation usually chooses as preferential orientation either the c axis 
(correctly described then by the Ising model) or its orthogonal plane (an XY model 
is then correct). 

The three-dimensional site-diluted Heisenberg model correctly describes the ex- 
perimental behaviour of a large number of real dilute magnetic materials, see Table HT 
so we will be able to compare our numerical results with some experimental estim- 
ates. 

In this model, according to the Harris criterion |23], see Appendix 13 the disorder 
is irrelevant. We want to check this point through numerical simulation by measuring 
critical exponents and different cumulants for different values of the dilution. If they 
do not depend on the dilution and agree with the pure case values, they will all belong 
to the same Universality Class (UC) and the Harris criterion will be re-verified. 

In addition we will study the self-averaging properties of the model computing 
at criticality the quantity R^, which will be defined below and is a measure of the 
self-averageness of the susceptibility. We will show results strongly supporting that 
this cumulant is zero at the critical point, but only taking into account the 
scaling corrections. This runs against some theoretical predictions [28] but supports 
others [27l|29]. 

We will obtain high-precision measurements of the observables for each lattice 
size near the critical point, so it will be necessary to take into account their finite-size 
effects in order to obtain asymptotic results. This implies estimating the correction 
to scaling exponents, whose leading term is denoted u, related to irrelevant operators 
in the Renormalization Group (RG) language. To this end, we will use the shift of 
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Ref. 


Material 


7 


/3 


6 


|ll6jlQQ/l 


FeioNiroBiigSi 


1.387(12) 


0.378(15) 


4.50(5) 


|ll(jjlQq/l 


FeisNisTBiigSi 


1.386(12) 


0.367(15) 


4.50(5) 


|llfj]lQq/l 


Fei6Ni64Bii9Si 


1.386(14) 


0.360(15) 


4.86(4) 


|117llll8]ioo. 


FeaoNieoPuBe 


1.386(10) 


0.367(10) 


4.77(5) 


|117llll8]ioa. 


Fe4oNi4oPi4B6 


1.385(10) 


0.364(5) 


4.79(5) 


|119jipq7 


FegiZrg 


1.383(4) 


0.366(4) 


4.75(5) 


inn] 1997 


FeggCoZrio 


1.385(5) 


0.368(6) 


4.80(4) 


inn] 1997 


FegsCoaZrio 


1.389(6) 


0.363(5) 


4.81(5) 


|119jipp7 


Fe84Co6Zrio 


1.386(6) 


0.370(5) 


4.84(5) 


|120jigpp 


Fei.ssMni.isSi 


1.543(20) 


0.408(60) 


4.74(7) 


IT5D]ippp 


Fei.5oMni.5oSi 


1.274(60) 


0.383(10) 


4.45(19) 


|l^ljl999 


MnCr1.gIno.1S4 


1.39(1) 


0.36(1) 


4.814(14) 


|l2ljiqqq 


MnCr1.8lno.2S4 


1.39(1) 


0.36(1) 


4.795(10) 


|122j9nnn 


Fe86Mn4Zrio 


1.381(12) 


0.361 




jl22j9000 


Fe82Mn8Zrio 


1.367(12) 


0.363 




jl23j9ooi 


Fe84Mn6Zrio 


1.37(3) 


0.359 


4.81(4) 


|123j9nni 


Fe74Mni6Zrio 


1.39(5) 


0.361 


4.86(3) 



Table 4.1: Experimentally-obtained critical exponents of materials which are ex- 
pected to be described by the three-dimensional site-diluted Heisenberg model with 
quenched disorder. Table from Ref. IH]- The results we obtain in this work are: 
7 = 1.398(6), /3 = 0.370(2), and 6 = 4.775(5). 



the crossing points both for the Binder cumulant and for the correlation length for 
lattice pairs of different sizes near the critical point. This study will also provide 
estimates of the asymptotic critical temperature value. We will also check that 
including the correction to scaling terms is crucial for the comparison of the values 
we obtain for the critical exponents with those of other workers. 

The simulations of this chapter were done mainly on the BIFI cluster. This 
consists of Xeon Dual Core 64-bit 3.40 GHz processors, with 2 GB of shared RAM. 
We used around fifty nodes for nine months making a total of around 17 years of 
computation time. 

4.2 Analytical Framework 

The self-averaging (SA) of the susceptibility is defined in terms of: 

R^^^^^^^W , (4.1) 
(An 
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with M. being the total magnetisation. The susceptibihty is self-averaging if i?^ — )• 
as L — )■ cxD. 

In Ref. [27J, the following picture was found: 

1. Away from the critical temperature: i?^ = 0. On the basis of the RG or 
using general statistical arguments, one can find that R^^ oc {i/LY in a finite 
geometry, L being the system size and ^ the correlation length which is finite 
for T ^Tc. Then i?-^ — )■ as L — )■ oo. This is called Strong SA. 

2. At the critical temperature, a RG analysis opens up two possible scenarios: 

• Models in which according to the Harris criterion the disorder is relevant 
(ctpure > 0): 0. The susceptibihty at the critical point is not self- 
averaging. In particular, Ref. [27] shows that under these conditions R^ 
is proportional to the fixed-point value of the coupling which induces the 
disorder in the Hamiltonian, which controls the new UC. This is called 
No SA. 

• Models in which according to the Harris criterion the disorder is not 
relevant (opure < 0): R^ = 0. The susceptibility at the critical point is 
self- averaging. In a finite geometry R^ scales as L°'^'^ — )■ 0, where a and 
i> are the critical exponents of the pure system, which are the same in 
the disordered one. This is called Weak SA. 

The observable R^ has been measured in other dilute models, for example in 
the four- dimensional dilute Ising model, see Ref. [51]. In this model a Mean Field 
computation and a numerical one found a non-zero value for R^ although the dilute 
model was shown to belong to the same UC as the pure model, contradicting the con- 
clusions of Ref. [27] . One can claim that the logarithms involved in the upper critical 
dimension make the numerical analysis difficult. In particular it was found analyt- 
ically in the mean field that R^ = 0.31024 and numerically that R^ G [0.15,0.32]. 
Because of the logarithms, it was impossible to make an infinite volume extrapola- 
tion for the numerical values of R^. Notice that in this model the only fixed point is 
the Gaussian one (all the values of the couplings are zero) and, following Ref. [27], 
R^ should be zero. 

In addition a two-loop field theory calculation done in Ref. p8] predicts a non- 
zero value for R^ for the dilute Heisenberg model (in which the disorder is irrelevant, 
«pure = —0.134, see Ref. |124] ). The two- loop field theoretical prediction for a in the 
pure case was ctpure > 0, so that apparently this work is consistent with the findings 
of Ref. [27] . The starting point in Ref. [28] was the mean field computation done in 
Ref. [51], modifying it to take into account the vector degrees of freedom, introducing 
the fluctuations using the Brezin-Zinn- Justin (BZJ) method, Ref. [125j . They found 
analytically R^ = 0.022688 for the vector channel and universally (independent of 
the dilution for all p < 1). It is important to remark that in the BZJ method one 
fixes from the beginning the temperature of the system to the infinite volume critical 
value, working in a finite geometry, so in order to compute R^ in this scheme the 
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following sequence of limits is used: 

R* = lim lim RJL,T) , (4.2) 

where i?* is the infinite volume extrapolation at criticality of Rx{L, T), and Tc is the 
infinite volume critical temperature of the system. The other possible limit sequence 
that can be computed is: 

lim lim RAL,T), (4.3) 

which is zero even when the disorder is relevant since R^ oc L~'^ as T ^ T^. 

Hence, in order to test these discrepancies we simulated numerically the site- 
diluted three-dimensional Heisenberg model computing i?* in the vector and tensor 
channels. To perform this programme, in particular in doing the infinite volume 
extrapolations of cumulants and exponents, a proper use of the corrections to scaling 
is really important. 



4.3 The Model 

The Heisenberg site-diluted model in three dimensions is defined in terms of 0(3) 
spin variables placed at the nodes of a cubic three-dimensional lattice, with Hamilto- 
nian 

i^ = -/3$^e,e,5i-5,-, (4.4) 

<i,i> 

where the Si are three-dimensional vectors of unit modulus, and the sum is extended 
only over nearest neighbours. The disorder is introduced by the random variables 
ej which take value unity with probability p and zero with probability 1 — p. An 
actual {ej} configuration will be called a sample. 

In addition, as done in Ref. |126] . we define a tensorial channel associated with 
the vector S through the traceless tensor 

rf^ = 5r5f-V, a,/3 = l,2,3. (4.5) 
We define the total nearest-neighbour energy as 

S = '^titjSi ■ Sj , (4.6) 

and the normalised magnetisation for both channels as 



M 



^5^6.5,, (4.7) 
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with V = and L is the hnear lattice size. Because of the finite probabil- 
ity of reaching every minimal value for the free energy, the thermal average of 
Eqs. (14. 7p and (14. 8 p is zero in a finite lattice. Therefore, we have to define the 
order parameters as the 0(3) invariant scalars 



^VA^) , Mt= UtiMiy (4.9) 



M 



Notice that the mean value of a non- invariant 0(3) observable is automatically zero. 
We also define the two susceptibilities as: 



X = V{M^) , XT = V{tTMl). (4.10) 
A very useful quantity is the Binder parameter, defined as 



,:.r-\m,_ . (4.11) 

Another kind of Binder parameter, meaningless for the pure system, can be 
defined as 



y {M^Y - {M^) T {tiMlY - (tr-M^) ..... 

92 — 2 ' 92 — 2 ' \^-^'^) 

{M^) (tr7W|) 

and these are the quantities we use to estimate the self-averaging properties of the 
susceptibility (i?^) in both channels. 

A very convenient definition of the correlation length in a finite lattice is, see 
Ref. EH], 



e = ( T^TITTTTT ) ' (4-13) 



X/F-I 
,4sin2(7r/L) 

where F is defined in terms of the Fourier transform of the magnetisation 



^{k) = ^J2e"'-^erSr (4.14) 



V 

r 



as 

V 



F = - (|^(27r/L, 0, 0) P + 1^(0, 27r/L, 0) P + |^(0, 0, 27r/L) P) . (4.15) 

The same definition is also valid in the tensorial case. This definition is very well 
behaved for the FSS method we have employed, see Ref. |126] . Finally, we measure 
the specific heat as 

C = V-^WT^Wf ■ (4.16) 
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4.4 Numerical Results 
4.4.1 Methods 

The lattice sizes L we have studied are 8, 12, 16, 24, 32, 48, 64, and, only in the pure 
model, L = 96. We have simulated five values of the dilution apart from the pure 
case, p = 1. These values are p = 0.97, 0.95, 0.9, 0.7, and 0.5. 

Between each measurement of the observable described in Sec. 14. 3[ firstly, we 
update the spin variables using a Metropolis method over 10% of the individuals 
spins, chosen at random, then we perform a number (increasing with L) of cluster 
updates using a Wolff method - see Ref. This is our elementary Monte Carlo step 
(EMCS). The number of clusters traced (or Wolff updates) between measurements 
was chosen to yield a good value of the self-correlation time, see Ref. [7] , in our case 
always 1 < r < 2 (r being the integrated autocorrelation time of the energy, see 
Appendix ICj) . 

In order to work in thermally equilibrated systems, we perform a great number 
of EMCSs before starting measurements. We start the simulation always from ran- 
dom (hot) distributions of the spin variables, although we have checked that the 
averages do not change if we begin from cold configurations (i.e., all spins pointing 
in the same direction). In particular, we took 4 x 10^ measurements for the pure 
model, discarding about 10^ of the first measurements for L = 8 and increasing this 
number with the lattice size. For every lattice size, we performed 2 x 10^ quenched 
disorder realizations in the dilute models (except for p = 0.97 and p = 0.95 with 
only 10^ realizations) taking 100 measurements per sample after equilibration, in 
accordance with Ballesteros et al. [SI] who demonstrated that the best approach to 
minimising the statistical error is to simulate a great number of samples with just 
a few measurements in each one. 

To measure the critical exponents, we use the so-called quotient method [70] . 
which allows great statistical accuracy, see Appendix [Bl Therefore, firstly we needed 
to estimate by successive simulations the /3 point where 

a2L,p,p) ^ aL,p,p) ^ ^^^^^^ 

2Zy L 

for each pair of lattices (L,2L). Then we used re- weighting techniques to fine- 
tune this condition. These re-weighting techniques are used to (3 extrapolate the 
observables and calculate their (3 derivatives, always before the sample averaging is 
performed. The equations used are, see Appendix [Dl 

{0){/3 + A/3) = {Oe^^^)/{e^^^) , (4.18) 

df,'{d) = MO) = {OS-{0){S)). (4.19) 

These extrapolations are biased. For instance, the expectation value of equation 
fl4.19|) . when the averages are calculated with measurements is 

(l-|-)a,(0>. (4.20, 
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Hence, we have to correct this bias, see again Appendix [Dl An example of the effect 
of this correction is found in Fig. 14.11 a major bias affects the uncorrected numerical 
data, and the importance of taking this effect into account is clear. In addition, it 
is clear that the recipe of Ref. ^SlJ is working perfectly for A^^ = 100, which is the 
number of measurements per sample we have taken in this work. Therefore, we are 
very confident that all the data presented in this work are not biased due to the 
re- weighting technique. 



0.024 




=^ 0.040 



0.030 



Figure 4.1: The g2 cumulant in both channels for L = 64, p = 0.9 with 1000 
samples, /3simuiation = 0.79112, re-weighted at /3 = 0.79082 as a function of 1/Nm, 
with Njji being the number of measurements in each sample. We report data with 
Nm = 50, 100, 500, 1000, 5000, and 10000. The data without the bias correction 
proposed in Ref. [51J are marked with triangles while the corrected ones are marked 
with circles. We also mark with the dotted lines the selection used in this work 
(which corresponds to Nm = 100). Notice the importance of the correction of the 
bias if one performs re-weighting with the data. 



Also, we tried to use the solution for the bias obtained in Ref. |127j . where each 
sample is split into four parts, but the results were poor. This was due to the 
small number of measurements we take in each sample (10^), which leads to large 
differences between the averages in each quarter. 

To compute errors in the averages we used a jack-knife method, see Appendix O 
We defined twenty jack-knife blocks for the pure model in a single sample and one 
block for each sample in the dilute {p < 1) models. 

The calculated observables and critical exponents sometimes present, instead of 
a stable value, a monotonically decreasing one. For t], there is found this type of 
evolution with increasing L, but it is clearly weaker than for z/, see Tables I4.5H4.12] 
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In these cases an infinite volume extrapolation is called for. If hyperscaling holds, 
we expect finite-volume scaling corrections proportional to L~'^. This issue will be 
addressed in the next subsection. 



4.4.2 The Scaling Exponent uj 

As will be seen in Tables 14.51 to 14.121 in Sec. I4.4.4[ there are evident finite volume 
effects, especially for the thermal exponents and the cumulants [g^ and g2). So we 
have to use the equation 



xq 

V 



V 



ocL"", (4.21) 

(L,2L) 



which is a consequence of the scale hypothesis first derived in Ref. |128j . Con- 
sequently, choosing a good value for w is a crucial question. 

Exact results and RG calculations tell us that the disorder, being irrelevant in 
this model, induces scaling corrections with an exponent a/z/ ~ —0.188 (in V) [8]. 
In addition to this new scaling correction one must have that of the pure model, 
which is related to the coupling of the (0^)^ term in the Ginzburg-Landau theory. 
This exponent is assumed to be 0.8 |129[I130] (for the pure model). Hence, the 
leading correction is the exponent induced by the disorder. We will try to check 
this scenario by computing the 'leading' correction to the scaling exponent from the 
numerical data. 

First of all, we tried to estimate lo just by considering it as another tunable 
parameter in Eq. fl4.2ip applied to some physical quantities. In these fits, as a 
first approximation, we disregarded the possible correlations between the data for 
different L values. The results are presented in Table 14.21 If we perform a weighted 
averaging with these results we obtain u = 1.07(9) for the pure model and u = 
0.92(9), 00 = 0.81(7), and oo = 0.88(4) for the dilute model with p = 0.9,0.7, and 
0.5 respectively, in very good agreement with the value of the scaling correction 
exponent of the pure model. However, we think this method is not very consistent 
because of the variability of the results from one quantity to another as seen in 
Table Ol 

Another approach, following Ref. |131] , is to study the crossing points of scaling 
functions (such as ^/L and measured in pairs of lattices with sizes L and sL. 
The deviation of these crossing points from the infinite volume critical coupling will 
behave as 

A/3{L, sL) = f3{L, sL) - f3,{oo) oc ^ 7 -^"^"^ • (4-22) 



With this method we need an additional estimate for the thermal exponent u. We 
used, following |124] . the value z/ = 0.7113(11) for the pure model (notice the really 
small error in z/, so that we will discard it in the following), which is also a valid 
value for the dilute models because of the validity of the Harris criterion, and as can 
be checked with the data below. We fixed s = 2. In this approach, we only use the 
crossing points in the vectorial channel because they are cleaner. 
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o 


l^p=1.0 


l^p=0.9 




l^p=0.5 




1.45(52) 











Vmv 


1.62(80) 













— 





1.2 (1.1) 


0.68(46) 


Vmt 


— 





— 


0.73(46) 




— 





— 


— 




2.30(61) 





— 


0.62(47) 




— 





— 


— 




2.12(52) 


1.76(60) 


1.09(40) 


1.34 (27) 




1.08(21) 


1.21(31) 


0.61(12) 


0.45(10) 


e/L 






1.55(76) 


1.64(17) 


9l 


0.85(14) 


2.00(61) 


1.21(15) 


1.19(13) 


9l 


1.06(14) 




1.35(33) 


1.41(42) 
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0.81(16) 


0.89(9) 


0.94(7) 


T 
92 






0.63(12) 


0.72(10) 


^weighted 


1.07(9) 


0.92(9) 


0.81(7) 


0.88(4) 



Table 4.2: u values from the L — )■ oo extrapolations of some quantities. The last 
row gives the weighted average of each column. We disregarded data with error bars 
larger than 100% of the values themselves. Those disregarded data are shown in the 
table as — . 



Extrapolating these crossing points using Eq.f l4.22]) . we can plot the minimum 
of the of the fit as a function of u obtaining the upper part of Fig. 14.21 and the 
whole of Fig 14.31 To carry out these extrapolations we must take into account that 
the measurements of the crossing points are correlated in pairs, so that we have to 
use the definition that includes the whole self-covariance matrix 

N N 

xl = 5^5Z(xi-fit)(cov-i)i,„(x„-fit), (4.23) 

1=1 m=l 

with N being the number of crossing points, that is to say, the number of simulated 
L values minus two; xi is the value obtained for the observable x (in our case the 
coupling) at the crossing point for Li and 2Li, and "fit" is the value fitted to the 
form of Eq. fl4.22p (or to another scaling form) for L;. In addition 

(COV);^^ = {XmXl) - {Xra){xi) (4.24) 

can also be defined in terms of jack-knife blocks, see Ref. [7], as 

/V - 1 

(cov),^ = J2{xjf - {xi)){x'-f - {xj) , (4.25) 
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Figure 4.2: Top: as a function of uj deduced from the fits to L — )■ oo, Eq. fl4.22p . 
for the crossing point of ^/L and for the (L, 2L) pair for the pure model. Also 
shown is the combined x^? whose minimum is marked with the dotted line. Bottom: 
Extrapolated /3c(oo) as a function of u. The point where the two observables give 
the same extrapolated value is marked with the dashed line. 



where A*";, is the number of jack-knife blocks, are block variables, where the 

first subindex runs over L values while the second one runs over jack-knife blocks, 
and {xi) is the average of all block variables given L = Li. 

Also, following Ref. |131] . we can do a combined fit in u of the crossing points 
of / L and by defining 

Xcombined = X^^/L + XgJ' ; (4.26) 

using Eq. fl4.23p to calculate each of the right-hand-side terms and searching for the 
minimum of Xcombined- obtain the error in u by searching for the point Ui 

at which Xcombined(^i) = Xcombincd (^min) + 1, SO that the crror is Au = l^min - 
The results for these combined fits are shown in the upper part of Fig. 14.21 and in 
the whole of Fig. 14.31 With this method we find the values 

CO = 0.96(15) , (4.27) 

for the pure model and 

CO = 2.29(70) , 0.84(17) , 0.64(13) , (4.28) 

for the dilute models with p = 0.9, 0.7, and 0.5 respectively, in agreement with 
the value obtained in the pure model |124[ll26[ll29tll30] . except in the p = 0.9 case 
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Figure 4.3: Top: ^ function of u for the dilute {p = 0.7) model. Also shown 

is the combined x^- Bottom: ^ function of u for the dilute (p = 0.5) model. 



for which the value is two standard deviations away from u = 0.8 |1291I130] . One 
possibility is that we are computing the leading correction to the scaling exponent 
but with a large error. Another possibility is that in the p = 0.9 model the coefficient 
of the leading correction to the scaling vanishes or is very small. This result and 
the change in the slope of the data for p < 1 with respect to the p = 1 ones, as 
can be seen in Table 14.31 constitute evidence for the possible improved action found 
for p = 0.9, see Ref. ^127j . Therefore the u exponent that we are measuring in this 
case could correspond to the third irrelevant operator, instead of the second one 
(remember that following RG the first one is a/u —0.188). 

In addition, as also was done in Ref. |131] . we were able to estimate the correct 
value for a; as that producing the same /3c(oo) value for both the crossings of ^/L 
and g4, as can be seen in the lower part of Fig. 14.21 marked with the dotted line at 
oj = 0.88. This approach only works for the pure model in which such a point is 
found. With another p value the (3c{oo) estimates from C,/L and do not cross each 
other. 

In conclusion, we have shown that our data (for both the pure and the dilute 
models) are fully compatible with the value u = 0.80(1) obtained previously both 
numerically and analytically for the pure model 0. In addition, since the error bars in 
u are really small (1% of error) we have discarded the uncertainty in u in the analysis 
presented in this work. Since the error bars in the extrapolated quantities are much 

"'^ Field theoretical approaches (both fixed dimension and e-expansion) provide very accurate 
values for w. 0.782(13) and 0.794(18) (respectively) [129) . Recent numerical simulations provide 
the values 0.775(13) and 0.799(13) [HD] and 0.64(13) and 0.71(15) [HH] . 
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larger than the uncertainty caused by the error bars in oo, we fixed uj = 0.80. The 
extrapolations obtained in the rest of the chapter are all obtained using this value. 

Finally, it is interesting to note that in the analysis presented in this subsection 
we have seen no traces of the leading correction to the scaling exponent even for 
the strongest dilution we have simulated, which should be a/z^ ^ —0.188. One can 
explain this fact by assuming that the amplitudes of this scaling correction exponent 
are really small, so that we are seeing only the next-to-leading scaling correction. 

4.4.3 Self- Averaging of the Susceptibility 

Having checked that the value u = 0.80 describes the corrections to the scaling for 
both the pure and the dilute models, we can try to extrapolate the values of g2 to 
infinite volume. 

Numerical results for g2 and in both channels are presented in Table 14.31 for 
both pure (only g^) and dilute models. 



p 


L 


92 


T 

99 

J z 


gl 


T 


1.0 


8 








0.62243(4) 


0.5216(1) 




12 








0.62172(5) 


0.5189(2) 




16 
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0.0553(6) 


0.1138(13) 


0.6138(2) 


0.5085(5) 
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0.0474(5) 


0.0980(11) 


0.6151(2) 


0.5095(4) 



Table 4.3: Cumulants for the 0(3) model. The first column is the spin density p. 
All the cumulants are calculated at the crossing points of ^/L for L and 2L. The 
averages were computed using 10^ samples (except in the p = 1 case). 
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Table 4.4: Cumulants for the 0(3) model with high p values (very soft dilution). In 
this case the cumulants are computed averaging 10^ samples. 



First of all, we will try to check the non-zero g2 scenario with the correction to the 
scaling exponent fixed to that obtained in the previous section. We found that it is 
possible, using the form of Eq. (I4.2ip (performing a combined fit) to extrapolate the 
values of to a value (depending only on the channel) which is independent of the 
dilution, and near the analytical prediction of reference |2B]- However, simulations 
at dilutions p = 0.95 and j9 = 0.97 do not follow the scaling found for j9 < 0.90 
(see Table l474|) . Hence, as a whole our numerical data do not support the scenario 
g2 7^ 0, see Figs. 14.41 and 14.51 for the two channels. Notice, see also Table S31 that 
all the values for these two lowest dilutions are smaller than the extrapolated point 
and they are decreasing (for both channels and taking into account the error bars). 

Secondly, we will check the 5'2 = scenario. To do this, we extrapolate g2 using 
the form proposed in Ref. [27] ((72 ^ L"/'") hut also including the term L~'^, i.e. we 
fit to: 

g2 = aL""^" + bL-'^ . (4.29) 

We obtain the fits shown in Figs. 14.61 and 14.71 for the two channels. The of these 
fits are really good. Hence, we have obtained strong evidence supporting this (72 = 
scenario. Notice that the introduction of the two scaling correction exponents has 
been of paramount importance for obtaining a very good fo^' ^ill the fits. The 
numerical data, for the simulated lattice size, do not follow the one-term dependence 
g2 oc L"/'^. 

4.4.4 Critical Exponents and Cumulants 

In this subsection we will check the consistency of the u exponent obtained in the 
text by means of the computation of critical exponents and cumulants. In addition, 
we will check whether or not these sets of exponents are universal by comparing 
different dilutions with the pure model. In this analysis we will use the data for 
p = 0.9,0.7, and 0.5. 
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Figure 4.4: Combined extrapolation to L — )■ oo for the (?2 cumulant of the vectorial 
susceptibility. Extrapolations were carried out by choosing a common value for the 
first term of Eq. (14.211) for all dilutions and by minimising the combined x^- We 
disregarded the data with L = 8 to obtain a good value for the x^- 
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Equation fIB.Sp applied to the quantities dj3^, dpg^^ M, and x yields respectively 
the critical exponents l + l/u, l/u, {D — 2 + rj)/2, and 2 — 1]. Their numerical results 
are given in Tables 145) and for the pure model, Tables H?71 and ITS) for the p = 0.9 
case, Tables K9\ and ICTD for p = 0.7, and Tables KTT\ and KW for p = 0.5. We also 
carried out combined extrapolations for all p values by fixing the same value of the 
extrapolated exponents for every p value. Some of these fits are shown in Figs. 14.81 
to 14. m and the compared results are presented in Tables 14.131 and 14.141 

The combined extrapolation of the Binder cumulant is given in Table 14.151 
The agreement of our results with those obtained in Refs. [124] (numerical for the 
pure model) and [28] (analytical) is really very good. We also obtain complete 
agreement with previous numerical estimates of the pure model critical exponents, 
see Ref. pi] . 

We obtain non-universal critical exponents and cumulants if instead of u = 0.8 
we use u = —ajv as the correction to scaling exponent. In addition, the dilution 
dependent exponents and cumulants are clearly different from the pure ones. Fur- 
thermore, this scenario does not change if we fit the data using both —ajv and 
w = 0.8. 
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0.987 


0.950 
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0.582 



Table 4.5: Magnetic exponents for the pure 0(3) model. The last three rows cor- 
respond to the L — )■ oo extrapolation (disregarding data with L = 8). 
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Table 4.6: Thermal critical exponents for the pure 0(3) model. In the last column 
we have disregarded data with L < 16. 
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Table 4.7: Magnetic exponents for the dilute 0(3) model with p = 0.9 . Extrapola- 
tions were carried out without disregarding data. 
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Table 4.8: Thermal exponents for the dilute 0(3) model with p = 0.9 . In the second 
and fourth columns we obtain poor results because the series are not monotonically 
decreasing. 
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Table 4.9: Magnetic exponents for the dilute 0(3) model with p = 0.7 . Extrapola- 
tions were carried out without disregarding data. 
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Table 4.10: Thermal exponents for the dilute 0(3) model with p = 0.7 . In the third 
column the fit was obtained disregarding data with L < 16. 





V 


Vt 


L 


X 


M 


XT 


Mt 


8 


0.0505(45) 


0.0461(48) 


1.3435(61) 


1.3431(62) 


12 


0.0448(39) 


0.0439(42) 


1.3684(54) 


1.3702(56) 


16 


0.0421(36) 


0.0417(39) 


1.3877(51) 


1.3896(52) 


24 


0.0396(32) 


0.0406(35) 


1.4033(46) 


1.4053(48) 


32 


0.0399(30) 


0.0414(32) 


1.4126(43) 


1.4152(45) 


L — )■ oo 


0.0346(60) 


0.0378(46) 


1.446(12) 


1.449(12) 


xVd-o.f 


2.225/2 


2.191/3 


0.119/1 


0.327/1 


C.L. 
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Table 4.11: Magnetic exponents for the dilute 0(3) model with p — 0.5. In the 
fourth and fifth columns we disregarded data with L < 16. 
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Table 4.12: Thermal exponents for the dilute 0(3) model with p = 0.5. In the 
second column we only used data with L > 8 while in the third and fifth columns 
we only used data with L > 12. 
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0.0389(10) 


1.4251(13) 


1.4251(14) 
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6.675/12 


5.104/15 
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0.878 
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0.518 
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Ref. [124] 


0.0378(6) 









Table 4.13: Combined extrapolation with all p values for the magnetic exponent rj 
compared with the results from Ref. |124j . The first three rows correspond to our 
L — )■ oo extrapolation. 
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Figure 4.8: Combined extrapolation to L — > oo for the 77 exponent deduced from 
the vectorial susceptibility (x^)- Extrapolations were carried out by choosing a 
common value for the first term of Eq. fl4.2ip for all dilutions, and by minimising 
the combined x^- 
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Figure 4.9: Combined extrapolation with all p values to L — > 00 for the 77 exponent 
deduced from the vectorial magnetisation (M^). 
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Figure 4.10: Combined extrapolation with all p values to L — > oo for the v exponent 
deduced from dpg\. 
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Table 4.14: Combined extrapolation with all p values for the thermal exponent v 
compared with the results from Ref. |124j . 
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Table 4.15: Combined extrapolation to L — )■ oo with all p values for the Binder 
cumulant (74 defined in Eq. (14.111) . compared with results from Refs. [2B| and [124] . 
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Figure 4.11: Combined extrapolation to L — )■ oo for the v exponent deduced from 

4.5 Conclusions 

We have studied the critical properties of the site-dilute Heisenberg model for differ- 
ent values of the dilution. Our main aims were both to re-verify the Harris criterion 
and to check the self-averaging properties of the susceptibility. 

We studied in great detail the corrections to the scaling in the model, finding 
that the numerical data follow the next-to-leading correction to the scaling exponent 
instead of the leading one. We obtained all the critical exponents and cumulants 
using this next-to-leading exponent. Also, the result of this analysis was found to 
be fully compatible with the RG predictions and the Harris criterion: our exponents 
and cumulants are compatible with those of the pure model and independent of the 
dilution to a high degree of precision. 

In addition, we showed that we obtain non- universal quantities if we assume ajv 
to be the main scaling correction even if we add the w correction to the scaling 
exponent, using two correction-to-scaling exponents in the analysis. 

Finally, we showed strong evidence for a zero g2 cumulant, in both the vector 
and the tensor channels, in the thermodynamic limit at criticality, contrasting with 
some analytical predictions [2H], but in agreement with others The introduction 
of scaling corrections in the analysis was crucial to obtain the (72 = scenario. In 
addition, simulations of samples with very soft dilution (p > 0.9) helped us to discard 
the (72 7^ scenario. 
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Chapter 5 

The Site-Diluted Ising Model in 
Four Dimensions 

5.1 Introduction 

One of the major achievements of statistical physics is the fundamental explanation 
of critical behaviour at continuous phase transitions through Wilson's Renormaliza- 
tion Group (RG) approach. While this has mostly provided a satisfying picture for 
over thirty years, certain types of phase transitions have resisted full treatment. Such 
stubborn cases, which have been the subject of conflicting proposals and analyses, 
include systems with in-built disorder. 

The Ising model with uncorrelated, quenched random-site or random-bond dis- 
order is a classic example of such systems and has been controversial in both two 
and four dimensions. In these dimensions, the leading exponent a which charac- 
terises the specific heat critical behaviour vanishes and no Harris prediction for the 
consequences of quenched disorder can be made [23] , see Appendix |X1 In the two- 
dimensional case, the controversy concerns the strong universality hypothesis which 
maintains that the leading critical exponents remain the same as in the pure case, 
and the weak universality hypothesis, which favours dilution-dependent leading crit- 
ical exponents (see |132] and references therein). 

Since D = 4 marks the upper critical dimensionality of the model, the leading 
critical exponents there must be given by mean field theory and there is no weak 
universality hypothesis. However, unusual corrections to scaling characterise this 
model, and the precise nature of these corrections has been debated. This debate 
motivates the work presented in this chapter: methods similar to those employed 
in |132] . namely a high-statistics Monte Carlo (MC) approach coupled with finite- 
size scaling (FSS), are used to advance our understanding of the four- dimensional 
version of the random-site Ising model (RSIM). 

While not directly experimentally accessable, the four- dimensional RSIM is of 
interest for the following reasons: (i) it is closely related to the experimentally im- 
portant dipolar Ising systems in three dimensions, (ii) it is an important testing 
ground for the widespread applicability of the RG, (iii) it presents unusual correc- 
tions to scaling, (iv) in high energy physics, the establishment of a non-trivial Higgs 
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sector |133] for the standard model requires a non-Gaussian fixed point and a new 
universality class which may, in principle, result from site dilution, and (v) it is the 
subject of at least five analytical papers which differ in the detail of the scaling 
behaviour at the phase transition. 

5.2 Analytical Framework 

5.2.1 Scaling in the RSIM in Four Dimensions 

The consensus in the literature is that the following structure characterises the 
scaling behaviour of the specific heat, the susceptibility, and the correlation length 
at the second-order phase transition in the RSIM in four dimensions (up to higher- 
order corrections to scaling terms) [TFl - ISTIfTM] : 



Coo(t) 



Xoo(t) ~ 



Here, the subscript indicates the size of the system, the reduced temperature t = 
(T — Tc)/Tc marks the distance of the temperature T from its critical value Tc, and 
A and -B > are constants. The correlation function at criticality decays as [IHtlSO] 

g^x) = x-^^'-^^^^llogxl^ (5.4) 

where x measures distance across the lattice, the dimensionality of which is D. The 
correlation length for a system of finite linear extent L also exhibits a logarithmic 
correction and is of the form 

a(t = 0)~L(logL)^ (5.5) 

The leading power-law behaviour is believed to be mean field because the fixed 
point is expected to be Gaussian and therefore 

« = 0, (3 = ^, 7 = 1, ^ = 3, iy=^, 77 = 0, A = ^. (5.6) 

Here, /3 and 6 are, in standard notation, the critical exponents for the magnetisation 
out of field and in field respectively while A is the gap exponent characterising 
the Yang-Lee edge. There is no dispute in the literature regarding these leading 
exponents, some of which will be re-verified in this chapter. Neither is there any 
dispute regarding the details of the unusual exponential correction terms in (15. ip - 
(15. 3p . However there are at least five different sets of predictions for the exponents 



A-5|t|-"exp I -2^/ — |log|t|| l|log|t| 



\i\ ^exp I \l — \log\t\\ ||log|t||^. 



Itr^exp I ljA|iog|t|| lllogltir. 



(5.1) 
(5.2) 
(5.3) 
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of the logarithmic terms, which differ from their counterparts in the pure model, 
and a principle aim of this work is to investigate these predictions numerically. 

Aharony used a two-loop RG analysis to derive the unusual exponential terms 
in dLlD-dSil), and also found 07] 

a = i, 7 = 0, z> = 0. (5.7) 

In |1H] , Shalaev pointed out that Aharony's results needed to be refined and, by de- 
termining the beta function to three loops, gave predictions for the specific heat and 
the susceptibility which differ from those in [U] in the slowly varying multiplicative 
logarithmic factors: 

d = 1.2368, 7 = -0.3684, ?} = 0.0094. (5.8) 

Jug studied the a = line of n-component spin models in {n,D) space where D is 
the system's dimensionality, and thereby worked out the logarithmic corrections for 
the D = 4 n- vector model For the case at hand (n = 1), he obtained 

a = 1/2, 7 = 1/212^ 0.0047. (5.9) 



In [Sn], Geldart and De'Bell confirmed that to obtain the correct powers of | log 
the beta function has to be calculated to three loops, but the results of |50] differ 
from those of 08] in the powers of the logarithms which appear in the specific heat 
and in the correlation function: 

1.2463, 7 ^-0.3684, r/ = ^ = 0.0047. (5.10) 



Finally Ballesteros et al. [51] extended and corrected Aharony's computation to give 
the correction exponents: 

« = ^, 7 = Y^^ 0.0094, z> = 0, q=^. (5.11) 

So the detailed analytic scaling predictions of at least five groups of workers 
clash, and a number of questions arise: (i) Is each set of predictions self-consistent? 
(ii) What is the full set of predictions (i.e., extended to all observables) originated 
by each set? (iii) Can a simulation approach provide numerical support for the shift 
in the correction terms from their counterparts in the pure model? (iv) And can 
such a computational approach lend support to one or other of these five different 
sets of analytic predictions? Here the scaling relations for logarithmic corrections 
developed in |135[I136] are used to answer (ii), and it is shown that the answers 
to questions (i) and (iii) and to some extent (iv) are affirmative. In particular, 
numerical support is presented for the broad scenarios presented in [17 tH9ll5T] . 



Modification of the self-consistent scaling theory for logarithmic corrections of |135[ 
1136] to incorporate the exponential terms leads to the following forms for the beha- 
viour of the magnetisation in the 4D RSIM: 



moo(t) = t'^exp |^-i^^|log|t||j|logt|^ (5.12) 
m^{h) = h^\loghf. (5.13) 
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The Lee- Yang edge, denoted by riy(t), is related to the locus of the Lee- Yang 
zeros along the imaginary /i-axis, see Appendix [Fl and marks the end of their dis- 
tribution. From Eq.(15) of [135] , we also write for its scaling in the paramagnetic 
phase 

rYL(t) -t^exp ^-^y^|log|t||j|logt|^. (5.14) 

Besides the scaling behaviour of the Yang-Lee edge, defined in Eq. fl5.14p . we 
also consider the density of zeros which, for an infinitely large system, we write as 
5'oo(^), where r parameterises their locus along the imaginary /i-axis (assuming the 
Lee- Yang theorem holds). In fact it is more convenient to consider the integrated, 
or cumulative, distribution function of zeros, defined as 

Goo{r,t) = / g^{s,t)ds . (5.15) 

J r-YLit) 

Following the approach outlined in |135] . its critical behaviour can be determined 
as 

Goo(r) ~ exp " |j) ' ^ogrr^'-'^^i , (5.16) 

where the exponential term drops out by using the mean-field values 7 = 1 and 
A = 3/2. 

The scaling relations for logarithmic corrections in this 4D model are 

a = Dq-Du, (5.17) 

2/3-7 = Dq-Du, (5.18) 

/3{6-l) = M-7, (5.19) 

V = 7-z>(2-r/), (5.20) 

A = P-j. (5.21) 

These scaling relations are now used to generate a complete scaling picture from 
the fragments available in the literature This complete picture is given in 

Table [5?H where the exponents of the logarithmic correction terms are listed. Values 
for the exponents in boldface come directly from the reference concerned and the 
remaining values are consequences of the scaling relations fl5.17l) - fl5.21l) . Each of 
the five papers [17H^ is self-consistent in that the exponents given within them 
do not violate logarithmic scaling relations. However, there are clear discrepancies 
between each of the five papers. 

The presence of the special exponential corrections has recently been verified by 
Hellmund and Janke in the case of the susceptibility [134] . These exponential terms 
mask the purely logarithmic corrections, so in order to detect and measure the latter 
one needs to cancel the former. Certain combinations of thermodynamic functions 



^ The relation (I5.17P is modified to read a = 1 + dq — dP when a — and when the impact 
angle of Fisher zeros onto the real axis is any value other than 7r/4, which is not expected to be 
the case in this 4D model [136] . 
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Log 


Pure model 


Aharony [47^ 


Shalaev gS] 


Jug m 




Geldart 


Ballesteros 


exp 


'5111137] 










& De'Bell fSU" 


et al fST] 


a 


1/3 


0.5 


1.237 


0.5 




1.246 


0.5 


/3 


1/3 


0.25 


0.434 


0.252 




0.439 


0.255 


7 


1/3 





-0.368 


0.005 




-0.368 


0.009 


6 


1/3 


0.167 


0.167 


0.170 




0.170 


0.173 


V 


1/6 





-0.189 






-0.187 





V 








0.009 






0.005 


0.009 


q 


1/4 


0.125 


0.120 






0.125 


0.125 


A 





0.25 


0.803 


0.248 




0.807 


0.245 



Table 5.1: Theoretical predictions for the exponents of the logarithmic corrections 
to scaling for the pure Ising model in four dimensions and for its random-site coun- 
terpart. The latter exponents are listed in boldface if they come directly from the 
cited literature. The remaining values are extended from those of the literature 
using the scaling relations f l5.17p -( [5.2ip . 



achieve this, but it turns out that FSS does this also. FSS therefore offers an ideal 
method to determine the exponents of the logarithmic corrections numerically [132] . 

5.2.2 Finite-Size Scaling 

Fixing the ratio of ^oo(i) in (15. 3p and ^l(O) in (15. 5p to x, one has 



r^'expl ^J^|log|t|| l|log|t|r = xL(logL)^ (5.22) 



Taking logarithms of both sides, one obtains 



log|t|| ^ -logL, (5.23) 



which re-inserted into (I5.22p gives 



t ~ L-t(logL)-'exp|i-^|ilogLl|l + 0^^1^ (5.24) 



~ i-^(losi)-»exp(^|logLj{l + o(^)}, (5.25) 

having used the mean-field value (15. 6p for the leading exponent u and the logarithmic 
scaling relation (I5.17p . If d = 1/2, this recovers a result in for the FSS of the 
pseudo-critical point. 

Inserting (15.250 into (15. 3p recovers (15. 5p . as it should. The FSS's of the remaining 
functions are determined by inserting (15.251) into (15.11) to (15.31) and (I5.12p to (15.141) . 
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One finds 

a \ / 6 



Cl{0) ^A- B'Lf exp (2 + ^ ^logi^j (logL)"+^«-'^) , (5.26) 

where B' cc B is a. positive constant |17H5T] . Inserting the mean-field values a = 0, 
V = 1/2, one obtains the simpler form 



Ci(0) ^A-B' exp (^-2^ ^ log L j (log Lf . (5.27) 

Similarly, the FSS for the susceptibility is 

Xl(0) ~L^|logL|^-^(^-«) = L2|logL|S (5.28) 

where 

C = 7-2(i>-g) = ^a + 7- (5-29) 
The FSS for the Yang-Lee edge is 

ri(L) ~ L-v|logL|^+v('^-« = L-3|logL|^ (5.30) 

where 

p = A + K^-q) = -\a-\l. (5.31) 



Each of these also has sub- leading scaling corrections of strength 0{1/ \J\ogL) times 
the leading behaviour. One notes, however, that the unusual exponential terms, 
which swamp the logarithmic corrections in the thermal scaling formulae (15. 2p and 
(15.141) . drop out of their FSS counterparts fl5.28p and (l5.3Up . These are therefore 
ideal quantities to study the logarithmic corrections. The theoretical analytical 
predictions of each of the five sources in the hterature are now used to construct 
five possible FSS scenarios for the specific heat, the susceptibility, and the Lee- Yang 
zeros. While Jug did not calculate the critical correlator or correlation length in 4D, 
the FSS picture corresponding to [lU] can still be constructed through the scaling 
relations for logarithmic corrections. The FSS scenarios are listed in Table 15.21 

The remainder of this chapter is concerned with Tables 15.11 and 15.21 The primary 
objective is to verify that the exponents for the logarithmic-correction terms in the 
RSIM are indeed different from those of the pure model. Once this is established, 
one would like to determine which of the five sets of analytical predictions are 
supported numerically. From Table [5^ it is clear that present-day numerics can not 
be sensitive enough to distinguish between all five scenarios for the susceptibility 
or individual zeros. However, there are clear differences between the predictions 
from imiiniEl] and from [l8l|50] for the specific heat (Table 15. ip , and it will turn 
out that the numerical data is indeed sensitive enough to favour the former over the 
latter. 
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Exponent 


Pure 
model 


Aharony 

m 


Shalaev 

m 


Jug 

[49] 


Geldart & 
De'Bell ^SOj 


Ballesteros 
et al [51j 


Susceptibility ( 
Lee- Yang zeros p 


1/2 
-1/4 


0.25 
-0.125 


0.25 
-0.125 


0.255 
-0.127 


0.255 
-0.127 


0.259 
-0.130 



Table 5.2: The exponents of the multiplicative logarithmic corrections to FSS for 
the magnetic susceptibility and for the Lee- Yang zeros coming from the literature 
and compared to their equivalents in the pure case. The FSS exponents are ( for 
the susceptibility and p for the Yang-Lee edge. 



5.3 The Model 

The partition function of the RSIM in a reduced magnetic field h is 

hY,eia, , (5.32) 

Wi} \ (id) « / 

where L denotes the linear extent of the lattice, the sum over configurations {cTj} 
is taken over Ising spins cTj G {±1}, denotes nearest neighbours, and ej are 

independent quenched random variables which take the value unity with probability 
p and zero with probability 1 — p. Below the percolation threshold {pc = 0.197 
in four dimensions), the phase transition is expected to disappear, while for every 
p > Pc there exists a critical (inverse) temperature Pdp) for each given dilution. 

In order to find the Lee- Yang zeros we define the energy, E, and the magnetisa- 
tion, M, of the system as 

-E = - ^ ^i^jCTiffj , M = ^ eiai , (5.33) 

(ij) i 

and 

Pl{/3 ■,M) = Y, Pl{E, M)exp{-/3E) , (5.34) 

E 

where the spectral density Pl{E, M) gives the relative weight of configurations with 
given values of E and M, the partition function in an imaginary field ih is therefore 

ZL{(3,h) = J2pL{P;M)exp{ihM) = Zi(/3, 0)(cos(/iM) + z sin(/iM)) , (5.35) 

M 

where the thermal average ((■ ■ ■ )) is a real measure, i.e., it is taken with Z{(3, h = 
0). Assuming the Lee- Yang theorem holds |138[I139] . since odd moments of the 
magnetisation vanish in the paramagnetic phase, the zeros for a given realization of 
the disorder are given by the values of h for which 

(cos(/iM)) = 0. (5.36) 
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In this way we obtain the zeros of the partition function for each value of p and L. 
Then we average over reahzations of the disorder (samples), and the resulting jth 
Lee- Yang zero is denoted by rj{L), the zero with j = 1 being the smallest. 

A robust method to determine the density of zeros, defined in Eq. f lS.lSp . from 
simulation data was presented in Ref. |140j . Defining the density of zeros for a finite 
system of size L along the singular line r > VYiit) as 

^?i(r) = L-^^(5[r-r,(L)], (5.37) 
j 

we can insert it into the cumulative density of zeros to obtain 

GL{r) = ^ gL{s)ds = for: r,(L) < r < r,+i(L) , (5.38) 

so that it is given at a zero by the average 

G'L[r,(L)] = ^. (5.39) 

We also measure the non-connected susceptibility, xwi defined as 

XL = ^{M^). (5.40) 

with V = being the volume of the system. This quantity is directly related to the 
average size of the clusters constructed using a Wolff algorithm [TD] . We checked this 
point in this work. In all cases, the two definitions of the non-connected susceptibility 
are fully compatible. 

Finally we measure the specific heat of the system, defined as 

C = ^{{E') - {Ef) . (5.41) 

5.4 Numerical Results 
5.4.1 Methods 

We performed extensive simulations of the model for linear lattice sizes from L = 8 
to L = 48 at dilutions p = 1, p = 0.8, and p = 0.5. In each case, we employed a Wolff 
single-cluster algorithm [TD] to update the spin variables using periodic boundary 
conditions. Thermalization tests including the comparison between cold (all spins 
up) and hot (all spins random) starts were carried out. We found that the plateau 
for the susceptibility is quickly reached by starting from cold configurations, see 
Fig. 15.11 Indeed, the results for the susceptibilities from hot and cold starts are fully 
compatible (and are less than two standard deviations away from each other, even 
at the level of logarithmic corrections). The information about the numerical details 
is given in Table 15.31 We took 1000 disorder realizations in all the cases except for 
L = 48, where only 800 samples were used. We estimate that the total simulation 
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Figure 5.1: Averaged behaviour of the susceptibihty with the MC time for 20 samples 
at L = 32 and p = 0.800. Measurements were performed after every MC sweep 
(Wolff update). The plateau is reached more easily starting from a cold configura- 
tion. 



time was equivalent to 20 years of a single node of a Pentium Intel Core2 Quad 
2.66 GHz processor. Since our aim is to estimate the scaling of quantities right 
at the critical point, simulations must be performed at the critical temperature of 
the model. We used the estimates for the critical temperature given in [SI]. In 
terms of /3 = 1/kT, where k is the Boltzmann constant, these are /3c = 0.149695, 
= 0.188864, and /3c = 0.317368, for p = 1, p = 0.8, and p = 0.5, respectively. 
In addition we simulated the dilution p = 0.650 at /3c = 0.235049 [51] using 
the same statistics as for the other dilutions. In this case we found the behaviour 
of the observables to differ from the expected. For example. Fig. 15.21 shows the 
strong deviation of the leading scaling behaviour of the susceptibility compared with 
that of the other dilutions. We re-checked this point starting from different initial 
configurations and even using different random number generators. This deviation 
is surely due to a biased estimate of the critical temperature in [SI]. For this reason 
we omit p = 0.650 from our analysis. 



5.4.2 The Pure Case p=l 

To establish confidence in the present approach, the pure system is analysed first 
to test whether the method employed successfully quantitatively identifies the log- 
arithmic corrections which are well established there. 
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Spin Concentration 


L 


^Wolff 




p= 1.000 


8 


200 


2 


(/3c = 0.149695) 


12 


400 


8 




16 


1600 


32 




24 


2000 


128 




32 


3000 


400 




48 


4000 


1600 


p= 0.800 


8 


100 


1 


(/3c = 0.188864) 


12 


200 


4 




16 


800 


16 




24 


1000 


64 




32 


1500 


200 




48 


2000 


1250 


p= 0.500 


8 


100 


2 


(^c = 0.317368) 


12 


200 


8 




16 


800 


32 




24 


1000 


128 




32 


1500 


512 




48 


2000 


1250 



Table 5.3: Simulation details for each spin concentration p and system size L. Here, 
^Woiflf denotes the number of Wolff updates between consecutive measurements, 
is the number of dropped measurements at the beginning of the MC history (in 
units of 10^). We always performed 10^ measurements within each sample after 
equilibration. 
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Figure 5.2: Scaling of the susceptibility comparing the p = 0.650 case at f3c = 
0.235049 ^Slj with the other dilute cases (also simulated at their corresponding (3c s 
obtained from [51]). There can be seen the strong deviation for p = 0.650 from the 
expected leading behaviour x ~ . The point size is in every case larger than the 
corresponding error bar. 



The scaling and FSS of the pure model {p 
The specific heat FSS behaviour is given by 



1) are well understood [ ^1137] . 



Ci(0)~(logL)°~(logL) 



1/3 



(5.42) 



up to additive corrections. Fitting to this form for a over the full data set 8 < 
L < 48, one finds the estimate a = 0.42(4) with a goodness of fit given by a 
X^/d.o.f. = 3.9/3, C.L.=27%. The estimate is two standard deviations away from 
the known value 1/3. As elsewhere in this analysis, inclusion of sub- leading scaling 
correction terms in the fits does not ameliorate this result, which is similar to that 
reported in [51]. 

The FSS for the susceptibility is given in fl5.28|) with ( = 1/2. Fitting to the 
leading form 

Xl(0) ~ Li (5.43) 

gives 7/z/ = 2.16(1) for 8 < L < 48 and 7/z/ = 2.13(2) for 12 < L < 32, the 
difference from the theoretical value 7/z/ = 2 being ascribable to the presence of the 
logarithmic correction term. Accepting this mean- field value for 7/z/ and fitting to 



XL(0)~L2(logL)' 



(5.44) 
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gives the estimate ( = 0.48 ± 0.02 in the range 8 < L < 48, albeit with xV^-O-f. = 
12.3/3, C.L.=1%. 

The FSS for the individual Lee- Yang zeros is given in fl5.30p with p = —1/4 in 
the pure case. Fitting to the leading form 

rj{L) ~ L-v (5.45) 

gives A/u = 3.074(5) for 8 < L < 48, the difference from the theoretical mean-field 
value A/u = 3 being due to the corrections. Accepting this value and fitting to 

rj{L) {log LY (5.46) 

gives p = —0.22(2) in the range 8 < L < 48. This estimate is compatible with the 
known value p = —1/4. As one would expect, the higher zeros yield less accurate 
estimates (as they are further from the real simulation points) with p = —0.18(3) , 
p = —0.17(7) , and p = —0.10(14) from the second, third and fourth zeros respect- 
ively. These estimates are listed in Table 15.41 

Having established that the numerics give reasonable agreement with the pure 
theory at the leading and the logarithmic levels, we now perform a similar analysis 
in the presence of disorder. 

5.4.3 The Dilute Cases p = 0.8 and p = 0.5 

Since the ansatz f l5.27p for the specific heat in disordered systems is somewhat more 
complex than that for the pure case f l5.42p . we begin the p 1 analyses with the 
susceptibility and the Lee- Yang zeros. It will turn out that our analyses will reinforce 
the analytical predictions that scaling is governed by the Gaussian fixed point and 
that the logarithmic corrections in the RSIM differ from those in the pure model. 
Indeed, the results for the zeros will be seen to be broadly compatible with the 
analytic predictions contained in |T7l-[5T]. 

For the weaker dilution value p = 0.8, a fit using all lattice sizes to the leading 
form fl5.43p for the susceptibility yields the estimate ^jv = 2.14 ± 0.01. Ascribing 
the difference from the Gaussian value 7/z/ = 2 as being due to the correction terms 
and, as in the pure case, fitting to fl5.44p . one finds an estimate for the correction 
exponent C = 0.39(3) for 8 < L < 48. This value lies between the pure value C, = 0.5 
and the theoretical estimates for the dilute value which are ( ~ 0.25 to 0.26. Thus, 
while the FSS for the susceptibility does not capture the theoretical estimates for 
the dilute case, the fitted values have moved away from the pure value and towards 
the lower value listed in Table 15.21 As elsewhere in this work, the inclusion of scaling 
corrections does not alter these results significantly. 

A similar analysis for the FSS of the susceptibility at the stronger dilution value 
p = 0.5 gives similar results: the leading form fl5.43p yields an estimate j/u = 
2.13 ± 0.02 with a goodness of fit given by x^/d.o.f. = 1.2/3, C.L.=75%. Ascribing 
the difference from the mean- field value 7/z/ = 2 as being due to the logarithmic 
corrections, and fitting to (15.441) . one obtains the estimate ( = 0.37(4) for 8 < L < 
48. Again this result is between the theoretical predictions for the pure {( = 0.5) 
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-0.18(3) 


-0.17(7) 
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0.8 


L = 8 - 48 


0.39(3) 
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-0.16(3) 


-0.20(3) 


-0.17(3) 


0.8 


L = 12 - 48 


0.42(4) 


-0.17(4) 


-0.16(4) 


-0.17(5) 


-0.18(4) 


0.5 


L = 8 - 48 


0.37(4) 


-0.20(4) 


-0.22(4) 


-0.21(4) 


-0.21(4) 


0.5 


L = 12 - 48 


0.40(6) 


-0.16(5) 


-0.20(5) 


-0.18(5) 


-0.19(5) 



Table 5.4: FSS estimates for the various dilution values, using a range of lattice 
sizes. For p < 1 the susceptibility is expected to scale as xl ~ -^^(logi^)^ and the 
Lee- Yang zeros as rj ~ L~^(logL)^, where ( ~ 0.25 to 0.259 and p ~ —0.125 to 
—0.130. (For comparison, the pure theory with p = 1 has ^ = 1/2 and p = —1/4.) 




2 2.5 3 3.5 4 0.8 1 1.2 1.4 

log L log (log L) 

Figure 5.3: Left: FSS plot for xl p = 0.8 (circles) and p = 0.5 (triangles) at 
the critical point. The slopes of the fitted solid and dashed lines give estimates for 
7/z/ of 2.14(1) and 2.13(2), respectively. Right: Plot of logXL — 21ogL against 
log (log L) at p = 0.8 (circles) and p = 0.5 (triangles) giving slopes 0.39(3) and 
0.37(4), respectively, indicating slow crossover of multiplicative logarithmic correc- 
tions from the pure case (where ( = 0.5) to the dilute case, where the theoretical 
value is C ~ 0.25. 



and the dilute {( ~ 0.25 to 0.26) cases. These results are summarised in Table l5^ 
together with the results for the same fits but with the smallest lattices removed. 
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Since in each of the dilute cases the susceptibihty results lie between what is 
expected for the pure and for the dilute theories, we appeal to the Lee- Yang zeros 
since they are expected to provide a cleaner signal. 




2 2.5 3 3.5 4 0.8 1 1.2 1.4 

log L log (log L) 

Figure 5.4: Left: FSS plot for the Yang-Lee edge at p = 0.8 (circles) and p = 
0.5 (triangles). The slopes of the fitted solid and dashed lines give estimates for 
A/u of 3.055(4) and 3.07(2), respectively. Right: Plot of logri -t- 31ogL against 
log (log L) at p = 0.8 (circles) and p = 0.5 (triangles). Fits in the range L = 12 to 
L = 48 (plotted) give slopes —0.17(4) and —0.16(5), compatible with the literature 
predictions that range from p ^ —0.125 to p ~ —0.13. (For comparison, in the pure 
model, p = —1/4.) 

The leading behaviour is first examined by fitting each of the first four Lee- 
Yang zeros to Eq. fl5.45p . For the weaker dilution given by p = 0.8, one obtains 
A/iy = 3.055(8), 3.056(9), 3.069(11), and 3.060(10) from fits to the first, second, 
third, and fourth zeros, respectively, using all lattice sizes. The equivalent results 
for the stronger dilution value p = 0.5 are A/u = 3.068(13), 3.071(15), 3.072(12), 
and 3.071(11), respectively. All fits are of good quality with acceptable values of 
X^/d.o.f., which we refrain from detailing. Again, these are interpreted as being sup- 
portive of the mean- field leading behaviour A/u = 3 with logarithmic corrections. 

The logarithmic-correction exponents are estimated by fitting to Eq. fl5.46p . with 
the various theories indicating that p = —0.125 to —0.13. The strongest evidence 
supporting this comes, as it should, from the first zero (the Yang-Lee edge) for 
p = 0.8, which yields the estimate p = -0.15(2) (with xVd-o.f. = 1.8/3, C.L.=61%). 
As in the pure case, and as expected, estimates for p deteriorate as higher-index 
zeros are used. Dropping the smallest lattices from the analysis, however, leads to 
these estimates for p becoming more compatible with |17H^ . These results are 
summarised in Table 15.41 

The equivalent analysis for the stronger dilution value p = 0.5 is less clear, with 
an estimate p = —0.20(4) coming from the first zero when all lattices are included in 
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the fit (with x^/d.o.f. = 3.3/3, C.L.=35%). Dropping the smallest lattices, however, 
gives p = —0.16(5) (with x^/d.o.f. = 1.8/2, C.L.=40%), closer to the values coming 
from [TtHST] . Similar results are obtained for the higher zeros, and these are also 
summarised in Table l57il 

As a final check of the reliability of our results, we used the spectral energy 
method [14211143] to re-weight the data obtained at f3c to f3c ± /^/3c (taken again 
from [51]), finding that the new data sets are fully supportive of the previous result^. 

Having established that the leading FSS behaviour corresponds to that origin- 
ating from the Gaussian fixed point, and that the logarithmic corrections to scaling 
are different from those in the pure model and moreover are (at least in the case of 
the Yang-Lee edge) broadly compatible with the literature predictions |I7H5T] . we 
now attempt to distinguish between these broad predictions. To this end we turn to 
the specific heat. 




10- 



10 20 ^ 30 40 50 

Figure 5.5: The specific heat for p = 0.8 (circles) and p = 0.5 (triangles). The error 
bars are in every case smaller than the point size. The solid and dashed curves are 
best fits to the ansatz Eq. fl5.27p . with a = 1/2. 



Having established confidence in the validity of the mean-field values 7 = 1 and 
Z\ = 3/2 for the 4D RSIM, we may use the scaling relation a = 2 — 2 A + 7 to also 
establish the mean-field value a = 0. The ansatz Eq. (15.271) for the specific heat 
may now be used. This contains information which can be used to discriminate 
between some of the scenarios in the literature. In Table 15. 1[ one observes that 

^ We followed the recipe given in Appendix |D] to perform the extrapolation to infinite number 
of measurements per sample. 
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there is a striking difference between the estimates for the specific heat logarithmic- 
correction exponent a coming from [ISKSO] and from [ITJIHIEI] ■ While the former 
have relatively large values of a, the latter agree on a = 0.5. The simulated values 
of the specific heat at p = 0.8 and p = 0.5 are plotted in Fig. 15.51 The slope of the 
full specific heat curve (I5.27P is 

This vanishes when Cl{0) = A and when ^/\ogL = dA/53/12. The first of these 
is the asymptote L ^ oo, from which A can be determined for each dilution. The 
second occurrence of zero slope is for quite small lattice sizes, i.e., beneath lattice 
size L = 8. Therefore a ^ ■^12/53\/log8 ~ 0.7, which excludes the values a ~ 
1.237 and a 1.246 given in [MlEn]- In fact, a best fit to the ansatz Eq. 
gives A = 89(34), B' = 142(71), and a = 0.57(14) for p = 0.8, and A = 60(25), 
B' = 102(54), and a = 0.76(3) for p = 0.5. Fixing d = 1/2 in each case gives 
A = 76(4), B' = 115(13) for p = 0.8, and A = 22(3), B' = 19(9) for p = 0.5. These 
curves are plotted along with the specific heat measurements in Fig. 15.51 However, 
fixing the correction exponent a to the value given in [ISlinn] yields a best-fit value 
of B' which is negative in each case, contradictory to llTJEn]- Thus we can deem 
these values to be unlikely. 



5.5 Conclusions 

Numerical measurements of the leading critical exponents in the 4D RSIM have 
been presented, confirming that the phase transition in this model is governed by 
the Gaussian fixed point. We then turned to the corrections to scaling, for which 
there exist five distinct sets of predictions in the literature |T7l - [5T] . The scaling 
relations for logarithmic corrections were used to complete these sets, and their 
counterparts for finite-size systems were given. 

The measured values of the susceptibility FSS correction exponent, (, for the site- 
diluted model lie between the known value for the pure model and the theoretical 
estimates coming from |17H5T] for the disordered system. While this result illustrates 
slow crossover of the susceptibility, the lowest lying Lee- Yang zeros give a cleaner 
signal. The measured value for their logarithmic correction exponents were indeed 
found to be compatible with the theories. 

To discriminate between the five theories, the detailed finite-size scaling beha- 
viour of the specific heat was also examined. The analysis was clearly in favour of 
the analytical predictions of |Tr t E9 i l5T] over those of |18||50] . This was contrary to ex- 
pectation since the former involve only two loops in the perturbative RG expansion, 
while the latter take the expansion to three loops in the beta function. 
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We have presented in this work highly accurate numerical simulations of various 
models of phase transitions in the presence of dilution. We checked the validity of 
some recent work, being able to outperform their statistical accuracy. 

Firstly, we checked the goodness of a recently proposed microcanonical simula- 
tion method [13] that computes entropy, rather than free energy, to derive all the 
thermodynamic information. The results, both for the pure four-state {Q = 4) 
Potts model in two dimensions and for the pure Ising model in three dimensions, 
were fully compatible with the most recent canonical simulations [HlllH]. To achieve 
this, we applied (for the first time, to the best of our knowledge) the Nightingale 
phenomenological renormalization ^\ (also called the Quotient Method |70j) within 
a microcanonical ensemble. In addition we obtained the critical exponents in the 
microcanonical ensemble, checking the validity of the Fisher renormalization [58] for 
models with a constraint in the internal energy. 

Once we had set up a correct microcanonical simulation method, we used it to 
study the (inherently complicated) strong first-order phase transition of the three- 
dimensional Potts model with Q > 3 states. We confirmed that dilution dramatically 
smooths the transition up to a tricritical point, with spin concentration p^, where 
it becomes of second order [HHIB] . We were able to claim that is definitely less 
than unity, having obtained pc = 0.954(3) and pc = 0.922(1) in the Q = A and Q = 8 
cases, respectively. We also obtained that within the first order region in the Q = 4 
case the latent heat is a self- averaging quantity while the surface tension is not. As 
a future development we will study the scaling of latent heat, surface tension, and 
critical temperature within the first-order region for the Q = 8 case. In this case we 
will have to deal with not fully equilibrated systems, so that we will have to quantify 
the effects of the lack of thermalization. 

Within the canonical ensemble, we studied the critical properties of the Heis- 
enberg dilute model in three dimensions for different values of the dilution. Using 
the next-to-leading scaling correction, we obtained results fully compatible with the 
Renormalization Group predictions and with the Harris criterion: our exponents 
and cumulants in the dilute cases were compatibles with those of the pure model 
and independent of the dilution. We also obtained strong evidence for a zero g2 
cumulant, in both the vector and the tensor channels, in the thermodynamic limit 
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at criticality, contrasting with some analytical predictions [2B], but in agreement 
with others The introduction of scaling corrections into the analysis was found 
to be crucial to obtain the g2 = scenario. 

We also studied the site-diluted version of the Ising model in four dimensions, 
confirming that the phase transition in this model is governed by the Gaussian fixed 
point. The logarithmic corrections to scaling were analysed to try to discriminate 
between five distinct sets of predictions |17H5T]. The measured values of the sus- 
ceptibility logarithmic correction exponent in the dilute case lie between the known 
value for the pure model and all the theoretical estimates for the disordered sys- 
tem, indicative of a slow crossover to the dilute universality class. We were able to 
discriminate between the five theories by a detailed study of the finite-size scaling 
behaviour of the specific heat. The analysis is clearly in favour of the analytical pre- 
dictions of [Tr t HO t lST] over those of [181150] . Further theoretical effort should be made 
in this field because the favoured scenerio stems from computation up to two loops 
in the perturbative RG expansion, while the rejected scenario involves expansion up 
to three loops in the beta function. 

The numerical results of this thesis were only made possible by the intensive 
use of important supercomputing facilities. We obtained from their resources more 
than the equivalent of 400 years of computation time of a single last generation Pen- 
tium 2.5 GHz CPU. Specifically, we used the clusters at the "Instituto de Biocom- 
putacion y Fisica de Sistemas Complejos" (BIFI) and the "Barcelona Supercom- 
puting Centre" (BSC). In addition, we exploited the volunteer computing platform 
IBERCIVIS, for which we are in debt with all its developers and volunteers. 
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En el presente trabajo hemos realizado simulaciones numericas de alta precision 
de varies modelos de transiciones de fase en presencia de desorden. Con dichas 
simulaciones hemos logrado comprobar la validez de los trabajos mas recientes y 
hemos mejorado la precision de sus resultados. 

En primer lugar, hemos comprobado la validez de un metodo de simulacion den- 
tro del colectivo microcanonico propuesto recientemente |T3j. Dicho metodo utiliza 
la entropia, en lugar de la energi'a, para obtener toda la informacion termodinamica 
del sistema. Los resultados obtenidos, tanto para el modelo de Potts puro con cuatro 
estados {Q = 4) en dos dimensiones como para el modelo de Ising puro en tres di- 
mensiones, son completamente compatibles con las simulaciones dentro del colectivo 
canonico mas recientes [S1I13]- Para lograrlo hemos aplicado por primera vez, que 
nosotros sepamos, la renormalizacion fenomenologica de Nightingale [69] (tambien 
Uamada Metodo de los Cocientes [70]) dentro del colectivo microcanonico. Ademas 
hemos obtenido los exponentes cn'ticos microcanonicos comprobando la validez de 
la renormalizacion de Fisher [58] para un modelo con una ligadura en la energia. 

Una vez que hemos asegurado la bondad de nuestro metodo de simulacion mi- 
crocanonico, lo hemos utilizado para el estudio de la (inherentemente complicada) 
transicion de primer orden fuerte que tiene lugar en el modelo de Potts tridimensio- 
nal con Q > 3 estados. Hemos comprobado que la dilucion suaviza drasticamente 
la transicion hasta llegar a un punto tricritico, con concentracion de espines pc, 
donde la transicion pasa a ser de segundo orden jil] - H6] . Podemos afirmar que 
es definitivamente mas pequefia que la unidad, habiendo obtenido Pc = 0.954(3) y 
Pc = 0.922(1) en los casos Q = 4 j Q = 8 respectivamente. Tambien hemos obtenido 
que dentro de la region de primer orden para el caso Q = 4 el calor latente es una 
magnitud autopromediante mientras que la tension superficial no lo es. Como una 
futura investigacion se pretende estudiar el escalado del calor latente, la tension su- 
perficial y la temperatura critica dentro de la region de primer orden para el caso con 
Q = S, para ello tendremos que tratar con estados no completamente equilibrados 
por lo que tendremos que cuantificar los efectos debidos a la falta de termalizacion. 

Dentro del colectivo canonico hemos estudiado las propiedades cn'ticas del mo- 
delo de Heisenberg diluido en tres dimensiones para diferentes valores de la dilucion. 
Usando hasta segundo orden en correcciones de escala hemos obtenido resultados 
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completamente compatibles con las predicciones del Grupo de Renormalizacion y 
con el criterio de Harris: los exponentes y cumulantes obtenidos son compatibles 
con los del modelo puro e independientes de la dilucion. Ademas hemos obtenido 
evidencias importantes de un cumulante g2 nulo, tanto en el canal vectorial como 
en el tensorial, en el limite termodinamico y en el punto critico. Este ultimo resul- 
tado contrasta con algunas predicciones anali'ticas [2H] y concuerda con otras [27] . 
La introduccion de correcciones al escalado en nuestro analisis ha sido fundamental 
para obtener el escenario con g2 = 0. 

Tambien hemos estudiado la version con dilucion por sitios del modelo de Ising 
en cuatro dimensiones, confirmando que la transicion de fase de este modelo esta go- 
bernada por el punto fijo Gaussiano. Las correcciones logaritmicas al escalado fueron 
analizadas para tratar de discriminar entre cinco conjuntos distintos de prediccio- 
nes p71 - [5T] . Los valores medidos del exponente de correcciones logaritmicas de la 
susceptibilidad en el caso diluido se sitiian entre el valor conocido del modelo puro 
y todos los valores teoricos estimados para el modelo diluido, esto es un signo de la 
existencia un fenomeno de paso (crossover) muy lento hacia la clase de universalidad 
del modelo diluido. Hemos logrado discriminar entre las cinco teon'as en confiicto 
haciendo un estudio detallado de comportamiento de escalado con el tamano del sis- 
tema del calor especifico. El analisis claramente favorece las predicciones anallticas 
propuestas en [ITIIIHIET] sobre las propuestas en [l8l[50]. Un esfuerzo teorico adicio- 
nal parece necesario ya que el escenario favorecido procede de un calculo hasta dos 
"loops" de la expansion perturbativa del Grupo de Renormalizacion mientras que 
el escenario descartado procede de una expansion hasta tres "loops" . 

Los resultados numericos de esta tesis doctoral solo han sido posibles debido al 
uso exhaustivo de importantes infraestructuras de supercomputacion. Hemos obte- 
nido el equivalente a mas de 400 anos de tiempo de computo de un linico procesador 
Pentium 2.5 GHz de ultima generacion. Especfficamente hemos usado los clusters del 
"Instituto de Biocomputacion y Fisica de Sistemas Complejos" (BIFI) y del "Bar- 
celona Supercomputing Centre" (BSC). Ademas hemos explotado la plataforma de 
computacion voluntaria IBERCIVIS, por lo que estamos en una deuda profunda con 
los desarrolladores y voluntarios. 
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Appendix A 

The Harris Criterion 



Given the fact that real systems are almost always impure, it is crucial to quantify 
to what extent, if any, disorder affects their critical behaviour. A criterion, the 
so-called Harris criterion, makes it possible to predict quantitatively the effect of 
disorder by using the critical exponents of the pure system only [23]. According to 
this criterion, the impurities change the critical behaviour only if the specific heat 
exponent a of the pure system is positive (the specific heat of the pure system is 
divergent). In the opposite case, a < (the specific heat is finite), the impurities 
appear to be irrelevant, i.e. their presence does not alter the critical behaviour. We 
will derive the criterion following the approach of Ref . |141] . 

Let us consider a system with quenched disorder. This can be for example the 
presence of impurities at random sites in a crystal lattice. In the pure case, this sys- 
tem undergoes a continuous phase transition at a temperature Tc. This critical tem- 
perature is expected to change in presence of disorder because it introduces spatial 
inhomogeneities in the coordination numbei0. The thermodynamics of second-order 
phase transitions is dominated by large-scale fluctuations. The dominant scale, or 
the correlation length, ^ ~ \t\~^ goes to infinite as t ^ 0, with t = {T — Tc) /Tc being 
the reduced temperature. 

The strength of the disorder (in our example, the impurity concentration) is 
denoted by p, with p = being the pure case. As Tc is approached the following 
change of scale length takes place. First the correlation length of the fluctuations 
becomes much larger than the lattice spacing, and the system 'forgets' about the 
lattice. The only relevant scale that remains in the system is ^(t). When we move 
closer to the critical point, ^ grows and becomes larger than the average distance 
between impurities, so that the effective concentration of impurities, measured with 
respect to the correlation length, becomes larger. It should be stressed that such a 
situation is reached for an arbitrary small initial concentration p. If p ^ 1 there 
is no reason for believing that the effect of impurities will be small. 

We will discuss, for the sake of simplicity, the Harris criterion using a particular 
model: the D-dimensional Ising-like system described in terms of the scalar field 

"'^This is by definition the number of interacting neiglibours of a given spin. In a mean field 
calculation, one obtains that the critical temperature of the Ising model is Tc — 2qJ /ks, where q 
is the coordination number, J the coupling, and ks the Boltzmann constant. 
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Ginzburg-Landau Hamiltonian, see for example [141] : 



H= I d^x 



l{V<P{x)r+^-{t-6t{x))^\x) + \g^\x) 



(A.l) 



where the quenched disorder is described by random fluctuations of the effective 
transition temperature 6t{x) whose probability distribution is taken to be symmetric 
and Gausssian: 



mx)r , (A.2) 



where po is the normalisation constant. For notational simplicity, the sign of 6t{x) 
in Eq. f lA.ip is defined such that positive fluctuations lead to locally ordered regions. 



Configurations of the fields (f){x) that correspond local minima in H satisfy the 
saddle-point equation: 

- A(p{x) + t(j){x) + g(p^{x) = 5t{x)(f){x) , (A.3) 

Such localised solutions exist in regions of space where t — 6t{x) assumes negative 
values. Clearly, the solutions of Eq. ( 1A.3|) depend on a particular configuration of 
the function 6t{x) being inhomogeneous. Let us estimate under which conditions 
the quenched fluctuations of the effective transition temperature are the dominant 
factor for the local minima field configurations. 

Let us consider a large region f]^ of linear size L ^ 1. The spatially average 
value of the function St{x) in this region can be defined as follows: 

St{f2L) = ^ [ dx6t{x). (A.4) 

Correspondingly, for the characteristic values of the temperature fluctuations (aver- 
aged over realizations) in this region we get: 



5tL = ^JSt^inL) = . (A.5) 

Then, according to Eq. (1A.3P the average value of the order parameter ^(i^^,) in this 
region can be estimated from the equation: 

t + g(f)'^ = 6t{QL) ■ (A.6) 

One can obtain that if the value of t is sufficiently small, i.e. if 

5t(l?L)>t, (A.7) 



then the solutions of Eq. (IA.6P are defined only by the value of the random temper- 
ature fluctuation 

1/2 



0(^2^) ^ ± . (A.8) 
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Now let us estimate up to which sizes of locally ordered regions this may occur. 
According to Eq. (lA.Sp the condition dt^^ t yields: 



L«^. (A.9) 

On the other hand, the estimation of the order parameter in terms of the saddle- 
point equation f lA.6p can be correct only at scales much larger than the correlation 
length ^ ~ t^". Thus one has a lower bound for L 

L> r^ (A.IO) 

Therefore, quenched temperature fluctuations are relevant only when 

or 

t2-"^<p. (A.12) 

According to the Josephson scaling relation, a = 2 — uD, see for example pp. Thus 
one recovers the Harris criterion: if the specific heat critical exponent of the pure 
system is positive, then in the critical interval: 

t<t, =pi/" (A.13) 

the disorder becomes relevant. This argument identifies 1/a as the crossover ex- 
ponent associated with randomness. In this case, the critical exponents of the dis- 
ordered systems differ from those of thepure one. In particular, the value of a for 
the disordered systems is never positive q 

On the other hand, if the exponent a = 2 — vD < 0, the condition flA.131) can not 
be satisfied near Tc (at t <C 1), and therefore in this case a weak disorder remains 
irrelevant in the critical region. 

In the marginal situation, i.e. a = 0, which is the case, for instance, in the four- 
dimensional Ising model (Chapter [S]) or in the two-dimensional Ising model [132j . it 
can be demonstrated |141j that although the specific heat exponent in the disordered 
model remains zero, the forms of the logarithmic singularities are affected by the 
disorder. 



rigorous argument for a < in the disordered case, applicable to many situations, is given 

in [Tug. 
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Finite Size Scaling and the 
Quotient Method 



When doing numerical simulations we are restricted to finite systems and therefore 
we will never obtain infinite specific heats or susceptibilities at the critical point. 
Nonetheless, there exist different methods to study the critical behaviour of a phys- 
ical system, working with a finite number of degrees of freedom. 

Probably the most popular approach is the use of Finite Size Scaling (FSS) 
techniques. They are based on the study of the evolution of observables with the 
system size in order to obtain information about the behaviour of the system at the 
Thermodynamic Limit (TL). FSS is based on the scaling hypothesis |128] . which 
states that the behaviour of the system is governed by the ratio L/^{oo, t), where L 
is a characteristic scale of the system (for example its linear dimension) and C,{oo,t) 
is the correlation length of the infinite system. If this ratio is large, the system is 
basically in its TL; if it is small, it will be in the FSS regime. 

One of the consequences of the above statement is that the evolution of the mean 
value of a given observable, O, with the system size will follow the functional form 



with t = (T — Tc)/Tc being the reduced temperature and fo is a smooth function 
depending on the observable. The leading correction term exponent, w, is minus the 
eigenvalue of the first irrelevant operator of the theory, in terms of RG language. In 
the following we assume that we are in the critical region, so that L. Then in 
the last term we can neglect compared with . The above equation is one 
of the multiple forms of the FSS ansatz, which can also be derived from a pure RG 
analysis, see for example [7] or [56] for detailed explanations. 

There exist more practical forms of the ansatz. The observable O diverges in the 
TL according to: 



(0(oo,t)) 



/o(L/e(oo,t)) + 0(r^i^ 




(0(oo,t)) = t-"o + -- - . 
For the correlation length, = u, so that we can make the change 

t-o = oc e(oo,t)"o/''L-"o/'^L"o/'^ = g{L/i{oo,t))L 





(B.3) 
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and Eq. fIB.ip can be rewritten as 

(0(L, t)) = L^o/'^[FoiL/aoo, t)) + L-GoiL/aoo, t)) + ■ ■ ■ ] , (B.4) 

where it can be shown that the correction term is also a function of L/^(oo,t). 
Given that ^(oo,t) oc t~'^, the scale variable can be written in terms of the reduced 
temperature to obtain an alternative form of the ansatz: 

{0{L,t)) = L^'oZ-lFoitL^) + L-^Go{tL^) + ■ ■ ■] • (B.5) 

Moreover, since we can use Eq. (lB.4p for the correlation length and Fq is smooth, 
we can invert it to obtain ^(oo,t) as a function of C,{L,t), and thus arrive at the 
most useful form of the ansatz: 

(0(L, t)) = L^o/u^p^^L/aL, t)) + L--Go{L/im t)) + ■ ■ ■] . (B.6) 

All the quantities in the above equation can be measured on a finite lattice, so that 
this will be our starting point for the explanation of the quotient method |70j . 

If we form the quotient, Qo, of a given observable, O, measured for two lattice 
sizes Li = L and L2 = sL, with s > 1 , at just such a temperature that = s, or 
explicitly 

asL,t) aL,t) 



sL L ' 

the result will be the elimination of the scaling function Fq in Eq. flB.6l) , 



(B.7) 



Qo|g,=. = 5^°/^ + 0(^-^). (B.8) 

Typically one chooses s = 2. The critical exponent xq can easily be derived from the 
above equation. The fact that there are strong statistical correlations between the 
quotients in Eqs. (IB.7P and fIB.Sp can be used via a jack-knife method to decrease the 
statistical errors in the numerical estimates of critical exponents, see Appendix [0 
or Ref. [7J. 

In the present work we used the quantities: 

X ^ Xo = u{2-r]), (B.9) 
M ^ Xo = ^(D-2-r/), (B.IO) 

^ Xo = iy + l, (B.ll) 
dis94 ^ Xo = l. (B.12) 

In addition, the crossing points of the correlation length, i.e. the temperatures 
where the condition of Eq. (]B.7|) is satisfied, provide an estimate of the critical 
temperature of the transition. By applying Eq. (]B.5|) to the correlation length, 
assuming that the scaling functions are smooth, we can obtain for the (inverse) 
temperature of the crossing the following behaviour: 

(3t'-(3c^l^^L---^ . (B.13) 
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The method can be also apphed in a microcanonical context if a vahd FSS ansatz 
is available. In this case the role of the reduced temperature t is played by e — Cc, 
where e is the energy density of the system and Cc is the energy density at the critical 
point, see Sec. 12.2.31 for more details. 

The quotient method can be improved to speed up convergence if logarithmic 
corrections are present. In particular, if a given quantity, O, behaves in the TL as 



where z can be either the reduced temperature t or e — Cc, the critical canonical or 
microcanonical exponent calculated using Eq. (IB.Sp must be corrected according to: 



where we use primes to label corrected exponents. 

If we have enough analytical information about the logarithmic term exponents 
we can apply the correction exactly. For example, for the two-dimensional four-state 
Potts model the values of the logarithmic correction exponents are known analyt- 
ically [l3l[86]. Thus we can calculate the corrections accurately. The susceptibility 
behaves as 




(B.14) 




(B.15) 



X~L^/^(logL)-V8 



(B.16) 



and we easily arrive at 




(B.17) 



For the correlation length it is known that 






and therefore its temperature derivative scales as 



(B.19) 



resulting in a canonical exponent correction of 




(B.20) 



For the microcanonical u exponent, i^m, we use that 



e~L-^/2(logL)-=^/\ 



(B.21) 



so that 



5ee~^'/^(l0gL)3/^ 



(B.22) 



and 




(B.23) 
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Data Analysis: Autocorrelations 
and Error Estimation 

The goal of this appendix is to provide a brief resume of the main ideas for the 
data analysis of the output of a dynamic Monte Carlo (MC) simulation. For a more 
detailed study see for example Refs. [7t lll5j . Our aim is to describe the modern 
techniques that avoid the usual error sources in this kind of numerical study. 

Given that the output of a dynamic MC simulation is a sequence of system 
configuration^, {0}o, {0}i, {0)2, • • • , {4'}n, we have to take two crucial issues into 
account: 

1. The initial bias: We have to start every simulation from a physically unrepres- 
entative configuration (usually "hot", all the spins in random configurations, 
or "cold" , every spin in the same state). The first configurations are thus not 
representative of the equilibrium distribution (the Boltzmann weight). There 
will be an initial transient regime which must be discarded to avoid a system- 
atic source of error. If we discard the initial data in estimating the mean 
value {0)i3 at the inverse temperature (3, then: 

1 ^ 

(o>.«o-^Eo., (C.l) 

where we have distinguished the true mean value (O) from its estimate O. 

2. Error estimates in equilibrium: The output of every MC simulation must be 
a confidence interval around the estimated mean value. The true mean value 
must lie within this interval at a reasonable level of confidence if the cor- 
rect procedures have been applied. Once equilibrium is reached, correlations 
between consecutive system configurations make the statistical error a factor 
2Tint,o larger than that of the corresponding independent sampling case, where 
Tint,o is the integrated correlation time of the observable O, see below. 

"'^As one does not need to store all the configurations, but only the values of a few functions of 
them (observables) , what one really has is a sequence of numbers Oo,Oi,02, ■ ■ ■ ,On, where O is 
a generic observable. 
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Both these issues are related to the same object, namely the autocorrelation 
function. 



The Autocorrelation Function 



By definition, the equilibrium autocorrelation function of the observable O at time 
t is [7j: 



Cooit) = {OsO,+t)p-{0)l 



Y X 



-PH{Y) 



YX 



-(3H{X) 



-0{X) 



(C.2) 
(C.3) 



where: 



0{Y) is the value of the observable O for the system configuration Y . 

\T]yx is the probability of reaching the configuration Y starting from the 
configuration X in t steps; i.e., it is a sum over all possible paths connecting 
X and y in t steps. 

exp(— is the Boltzmann weight of the configuration F, with /3 being 
the inverse temperature, H the Hamiltonian, and Z the partition function. 



A normalised form is often used, defined as: 

Coo{t) 



pooit) = 



Coo(O) 



(C.4) 



Typically, for long times, Coo{t) decays exponentially with time. 
The exponential autocorrelation time is defined by 



T-exp,o = lim sup 



log Poo (i) ' 

It is useful to define the maximum over all the measured observables, 



Texp = sup rexp,0 
O 



(C.5) 



(C.6) 



In Ref. [115] it is demonstrated that the rate of convergence to equilibrium from an 
initial non-equilibrium distribution can be bounded in terms of Texp. In particular: 



J20{Ymyx-{0), 



~ e 



(C.7) 



From this, it can be said that setting = 20rcxp in Eq. fIC.ip is enough for all 
practical purposes. Hence, waiting for this time before starting to save the meas- 
urements for averaging, we can avoid the initial bias of the MC simulation. The 
problem with this approach is the difficulty in estimating Texp in some cases. 
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Usually the convergence to equilibrium is determined empirically by plotting 
certain observables as a function of time and noting when the initial transient seems 
to end. This includes the comparison between hot and cold starts. The main 
objection to this is the possibility of metastability, especially for first-order phase 
transitions. In such cases the equilibrium appears to be reached but really the system 
has just settled down into a long-lived metastable region of the configuration space. 
One has to be extremely careful in these cases, see Chapter [31 

Once in equilibrium, to what extent are the measurements taken in the system 
representative? This issue reflects the fact that consecutive measurements are usu- 
ally close in configuration space (and are thus correlated) so they do not provide the 
same information as if they were independent. 

We can resolve this question in terms of the integrated autocorrelation time, 
defined as: 



:^ + Y.Poo{t). (C.8) 



^int,0 

Z 

t=l 

This time controls the error estimates in MC simulations. In particular, the sample 
mean 

1 ^ 
t=i 

assuming for brevity that the data at t = 1 are already good, has a variance 

1 ^ 

{{0-{0),r), = ^5^((0.-(0),)(0.-(0),)), (C.IO) 

r,s=l 



1 ^ 

^2j2^ooir-s) (C.ll) 



iV2 

r,s=l 



1 / ifl\ 



t=-(N-l) 

{{0'),-{0)l). (C.13) 



27int ,0 f I r\2\ //n\2 



To derive these last relationships, we made use of Coo(^) = C'oo(— ^) in equilibrium 
and assumed that is large enough to neglect Coo(^) for \t\ ~ N - recall that 
Coo(^) exponentially for increasing t. 

Therefore the variance of O is a factor 2rint,o larger that it would be if all the 
configurations {0}j were statistically independent. In other words, the number of 
effective measurements in a MC simulation of length is reduced to iV/(2rint,o)- 
Roughly speaking, the error bars will be of order (rint,o/^)^^^5 so that if we want 
an error bar for our simulation of around 1% precision we will need a run of length 
^ lOOOOrint. 

Now we will define a more practical estimate of the correlation times |115j . The 
direct estimate from a run of length iV (supposing again that at t = 1 the data are 



139 
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already good) is: 



N-\t 



Coo{t) 



N-\t 



O) 



(C.14) 



Poo{t) = Cooit)/Coo{^)- (C.15) 
At first sight, one would estimate the integrated autocorrelation time as 



int,0 



1 ^-1 



(C.16) 



t=i 



But this is wrong because this estimator has a variance that does not go to zero for 
large A^. Each of the terms poo{t) decreases exponentially with t but is a random 
variable obtained by averaging N — t data. For t > Tint,o each pooif) is very small, 
but its error is not null. Therefore most of the terms in the definition Eq. (1C.16I) 
carry little information but much noise for large t, and there are very many of them. 
The cure is the truncation of Eq. (]C.16|) to estimate Ti^t,o self-consistently: 



int,0 



(C.17) 



t=i 



where a is a small fixed number, around 6. Clearly rint.o is biased, but the con- 
tribution of the neglected terms should be negligible (at 6rint,o, it is expected that 
Pooif) ^ !)• See |115] for more details about this "automatic windowing" al- 
gorithm. 
For r, 



exp,o, one can use two different estimates, see [7j: 

t 



log pooit)\ 



'exp,Ot 



log 



poo{t) 
Poo{t+l) 



-1 



(C.18) 



(C.19) 



The autocorrelation functions decay as pure exponentials, e~*/^'=''P'° , for large t. Then 
both rexp,Ot A ^"^^ Texp,Ot B become equal to rcxp,o- Nevertheless as was found before, 
the information carried by the signal decreases rapidly with t, and it can not be said 
which definition will first reach the t-independent region, see [7]. 



Error Estimation 

Having demonstrated that the effective number of measurements is A^/(2rint o)) 
where rint,o is not known accurately, we will now describe the most usual method of 
obtaining independent quantities: making data blocks. This will also allow the es- 
timate of functions of observables, and we will see that the time correlations between 
different observables can even be beneficial. 
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Starting from a set of data, let us form N/h data blocks of size h: 



^ bi 

Ob, = -^ Of (C.20) 

t={i-i)b+i 

The autocorrelation times for the blocked data are divided by a factor b. The mean 
value of the block data, O, is 6- invariant, while the error can be estimated (if b is 
large enough) as for statistically independent data: 







b ^ 



N/b N/b 
i=l \ i=l 



(C.21) 



The error estimate first grows with b until the data blocks become effectively inde- 
pendent (6 ^ 20rint o)- Then Aq becomes 6- independent (the expression in brackets 
in Eq (IC.2ip decreases as 1/6 but the prefactor grows as b). As can be seen in Ref. [7], 
an assumption that the data are independent would underestimate the errors by a 
factor ^^/2Tint,o- Once b is large enough, the error estimate reaches a plateau, but 
from then on the fluctuations grow as the block size is increased. Therefore N/b 
should be kept large enough. Typically, a ratio of between 10 and 50 is chosen. 

Now we have to consider the crucial issue of the error estimation of functions of 
observables. Let f{{0^^^)i3, {0^'^^)i3, . . . , (O^^^^)/?) be a function of R observables (R 
could be one). The best thing we can do is estimate /({(O*^^-*)/?}) by /({0(-^)}), 
although we know that this estimate is biased unless f{x) is a linear function. The 
estimate of the error could be done by linear error propagation, but when / depends 
on several observables this may overestimate the size of the errors due to correlations 
between the observables. For some clear examples of this last point, see Ref. [7]. 

The fact that correlations can be beneficial is exploited by the jack-knife method. 
This allows one to estimate the error bars of derivative quantities easily and coher- 
ently. The procedure is the following [7]: 



1. Estimate f{{0^'^)p, (O^^))^, . . . , {O^^))^) by f{OW, 0(2), . . . , 0(^)) . 

2. For each observable, form the corresponding block data, as in Eq. (]C.20p . with 
large enough b (the same for every observable). 

3. Make jack-knife blocks from the blocked data. This means that the i-th jack- 
knife block is formed by averaging all the blocks formed in the previous step 
except the i-th. I.e.: 



' r = l,2,...,R. (C.22) 

b j^i 

4. Estimate the function value for the jack-knife blocked observables as 

fjK,b,i = fiOjlfi^i, Ojl i^ i, . . . , O^r^ .) . (C.23) 
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Since the expression in brackets is an average of blocked data, it is smaller 
than usual. This is the reason for the multiplication (instead of division) by 
the number of blocks minus one. In the case of correlations between observ- 
ables, their jack-knife blocks will fluctuate simultaneously, thus reproducing 
the possible positive effect on the error. 



5. Estimate the variance of the function as 




(C.24) 
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Appendix D 

Temperature Extrapolations 



Within a canonical Monte Carlo method, the temperature of the system is kept 
fixed, and all the information about the observable corresponds to the simulation 
temperature. It is often very desirable to obtain an accurate estimate of a given 
quantity at a temperature different from the simulated one. This may be the case 
for example when one tries to obtain the absolute maxima in temperature of some 
quantities to estimate critical exponents, or when one has to fine-tune some condition 
as in the quotient method, see Eq. (lB.7p . 

Using the energy histogram of the system at a given temperature, one can ob- 
tain accurate information at nearby temperatures. The method was first proposed 
in |142j . and was recovered in |143] . If we are working with disordered systems, the 
temperature extrapolation must be performed before averaging over the different 
disorder realizations. 

The following formulae allow one to calculate the thermal derivative of an ob- 
servable, O, and its value at a temperature close to the simulated one: 



dp{0)=d^{0) = {OE-{0){E)), (D.l) 

where /3 = 1/T is the inverse simulation temperature. 

Nevertheless, it must be borne in mind that the two above expressions involve a 
systematic bias whose correction can become critical. We shall follow the approach 
proposed in jl44j . Using Eq. (]D.1|) with different measurements, we are really 
obtaining 



{^-^)mO), (D.3) 

where r represents the integrated autocorrelation time, see Appendix [Cl between 
the observable O and the energy. Since the r value is different for different samples, 
we shall have to obtain the disorder average. 

Let us demonstrate the validity of Eq. flD.3p . Consider two observables, and 
O^, and, using N measurements, calculate their connected correlation 

{O'^O'') ~ {0''){0'') , (D.4) 
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Since the mean value of O'^ is by definition 

1 ^ 

1=1 

Eq. f lD.4p can be written as 

^ N ^ N N 

NT.O^O\--,Y.Y.OtO], (D.6) 

where the latter term is more complex given that the measurements of the two 
observable may be correlated. We can rewrite it as 

N N N N-i 

J2 ot (^S = E E ou ■ (D-7) 

i=l j=l i=l t=l—i 

Recall that the behaviour of the correlation between two observables is for large 

\t\: 

(OtOl,) ^C, = {0-){0') + Coexp [-^ 



(D.8) 



where 

Co = (0^0^) - (0")(0^). (D.9) 

As we are summing over i, we can replace 0^_^_^ by its mean value, and in order 
to perform the sum we shall take 3> From the exponential decay with |t| 
in Eq. (ID.Sp we can extend the limit of the sum to infinity, and simplify things by 
replacing the sum with an integral. Therefore 



^ Of = {{O") (O'')) + NCo / dt exp - ' ' 

J — oo ^ 



i=i j=i 



(D.IO) 



But the integral is 2r, see again Appendix [Ul and replacing Cq by the expres- 
sion ( 1D.9P one recovers Eq. (ID.Sp . 

Thus it has been shown that the derivative of the mean value of an observable is 
subject to a systematic bias of order 2t/N, although corrections of higher order can 
be expected. This systematic error is not considered in general in other MC studies 
mainly because it is usually masked by the statistical error, of order I/a/ZV, that is 
far larger than the bias. 

Nevertheless, in a dilute model the situation is more complicated because there 
are two parameters involved: the number of measurements that we perform within 
each realization of the disorder (Ising average), denoted by Nj, and the number of 
realizations of the disorder (sample average), Nm- 

We can associate with every observable ctm, representing the deviation of this 
observable between different samples, and cx/, representing the average of the de- 
viations within each sample. Assuming statistical independence between different 
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samples (this is the case if there are no correlations in the random numbers) and 
independence between different measurements within a sample (possible using a 
cluster method), the variance of the mean value of an observable can be written as 



From this equation the optimal choice of Nj and Nm can be deduced if we want to 
optimise the simulation time to obtain a given error bar. Note that the simulation 
time is proportional to NmNj, because the simulation spends an approximately 
fixed portion of time taking measurements. For this reason the optimal value of Nj 
can not be much larger than a'j/a\j, i.e., Ni is limited by this latter expression. 
Therefore the best procedure to improve the statistical errors is to increase the 
number of samples, performing a not too large number of measurements within 
each one. Numerous studies on dilute models have avoided tackling this question 
by working in general in the regime Ni 3> Nm- Besides, a value of Ni that is large 
does not improve the result for the deviation a\j because it is only suppressed by a 
factor Nm- 

However, if we need to obtain the temperature derivatives there is a term in the 
statistical error proportional to 1/ ^/Nm- If this number is of order 1/Ni (systematic 
error in the extrapolation) a detectable bias appears. Unfortunately this is indeed 
the scenario in some cases. For example for the Heisenberg model, see Sec. 14.4. ![ 
we used Nm = 10^ and Ni = 100 measurements per sample. 

As a result it is necessary to find some algorithm to obtain correct results from 
the simulation data. The first possibility is to use fully independent measurements 
(by assuring 2r = 1) for each disorder realization. In such a case, recalling Eq. (ID.Sp . 
we could build a correct unbiased estimator by multiplying by the factor 1/(1 — l/Nj) 
in the case of the derivative and by the corresponding quantity in the case of the 
temperature extrapolations. However this is too expensive in simulation effort, due 
both to the need for too wide a time interval between measurements and to the 
fluctuations of r between different samples, which forces one to disregard too many 
measurements in samples with r smaller than the maximum. This makes the use of 
this approach inefficient. 

Another possibility, see Ref. [51], is to correct the systematic bias by splitting 
the measurements into statistically independent groups. In this way, by multiplying 
estimators from the different groups, correct values of {0°'){0^) can be obtained. 
Nevertheless, in a real situation independence is complicated because, in general, 
consecutive measurements in MC time are used to form these groups. 

In this work the approach that will be used for the dilute models, see again 
Ref. [51], is to extrapolate to A?^/ — )■ oo by using estimators coming from different 
values of Ni. We will work with 2r > 1, but the equality is not essential, as we will 
demonstrate in the following. 

In particular, for each disorder realization we can calculate the derivative with the 
entire MC history [Nj measurements), obtaining the estimator, whose systematic 
error is proportional to 2t/Ni. If we consider the two contiguous halves of the total 
history, we obtain two estimates for the derivative that, once averaged, produce 




(D.ll) 
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the estimator, y2, whose systematic error is proportional to 4t/Nj. If the history 
is divided into four, the derivative in each part is obtained and the average of the 
four parts is calculated, the estimator, y^, will have a bias of St/Nj. In Fig. ID.lt 
obtained in Ref. [Slj, is shown the observable dpx fo^' ^ four- dimensional site-diluted 
Ising model in a system with L = 48, p = 0.5. The estimator yi is plotted, but using 
only one of every k measurements of the history and averaging the results. It is 
obvious that doing so reduces the r by a factor k, as also is the case for the number 
of measurements, Nj. Therefore the bias, proportional to 2r/iV/, remains constant. 
This effect is maintained up to ~ 10. The correlation time between energy and %, 
r, was about eight measurements. 

Using the yi values, linear and quadratic extrapolations to Nj ^ oo in the 
variable 1/A^/ can be obtained with the form 

y. = y^ + ^^+^,. (D.12) 

For the linear case {B = 0) we obtain 

Vl = Voo = '2yi - y2 , (D.13) 

and for the quadratic case 

8 1 

yQ = yoo= gl/i - 2y2 + -ys . (D-14) 

The procedure is repeated for each sample and averaged over the disorder. The 
value yq is used to verify that the higher-order effects are negligible compared with 
the statistical error of the extrapolation in temperature. If this were not the case, 
to obtain measurements of the observable at the shifted temperature it would be 
advisable to make another simulation closer to that point. 

The estimate for the temperature extrapolations of the 0(3) model was obtained 
in this work in a similar way, see Chap. HI Fig. 14.11 
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Figure D.l: Upper part: Sample average of the estimator yi corresponding to dpx 
considering just one of each k measurements, for a four- dimensional Ising model on 
an L = 48 lattice, with p = 0.5 and t ^ 8. Lower part: Sample average of yi, y2, 
and 1/3 as a function of the inverse of Nj, on a lattice, with L = 32, p = 0.5 , and 
0.8. Figure taken from Ref. [51j . 
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Appendix E 

The Maxwell Construction 



Let us consider a system with action S[ip], depending on the local variables ip, 
coupled with a source J. If J) is the Gibbs free energ 30 and 0(x) = 0(^(x)) 
is a generic observable then 



where V is the system volume and D is the spatial dimensionality. The equation of 
state is 



relating the observable, O, with its conjugate variable, J. 

In the neighbourhood of the phase transition, the function O = 0{J) may be 
discontinuous. In other words, there may exist a range of O values that do not 
correspond to any J value, see Fig. lE.ll 

An effective potential associated with the observable O can be defined as: 






(E.2) 



e-vr(o) 



/ 



d[ip]e 




(E.3) 



where S is the Dirac delta function. 



Then it can be obtained 




(E.4) 



But for a large volume this integral 



is dominated by the saddle point 



d£ 
W 



J 



(E.5) 



and we can conclude that 0(/3, J) is the Legendre transform of F: 



^{f3,J) = r{0) + J0. 



(E.6) 



We have included in its definition the j3 term. 
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o 



Figure E.l: Lack of continuity of the function O = 0{J) at a first-order phase 
transition. 



In Eq. f lE.6|) . the O must be chosen that produces a global minimum of r{0) + 
JO. The contribution of all other minima will be exponentially suppressed. In 
general r{0) + JO has two local minima. One will be the absolute minimum for 
J < Jc, and the other will be that corresponding to J > Jc- Recall that O = 0{J) 
may be discontinuous. 

The condition for a minimum located at O is 

±{r{O) + JO)=0. (E.7) 

Let us define the two local minima as: 

min[r(0) + JO] = Oa , (E.8) 

min [r(0) + JO] = Ob ■ (E.9) 

J^Jc 



Therefore 

dr 



dr_ 
W 



-J,. (E.IO) 



Of 



Just at Jc both minima will result in the same value of r{0) + JO, i.e., 

r{OA) + JcO~A = r(OB) + JcO^. (e.ii) 

Given that _ 

r(0) = - [ dOJiO), (E.12) 

Jcte 
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Eq. flE.lip can be written as 



Ob 



dOJ{O)-.UOB~OA)=0 



(E.13) 



Oa 



which imphes that the shaded areas in Fig. IE.2I must be equal (in absolute value). 
This represents the well-known form of the Maxwell construction. 




Figure E.2: Usual form of the function J 
construction. 



J{0). Also shown is the Maxwell 
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Appendix F 
Lee- Yang Zeros 



In 1952 T. D. Lee and C. N. Yang, as part of their study of the phenomenon of 
spontaneous symmetry breaking, wrote two impressive papers |138] . In analysing 
the behaviour of a lattice gas (which is equivalent to an Ising model in a magnetic 
field), they approached the problem of its phase transition in an absolutely novel 
way focusing on finding the zeros of the partition function in terms of an external 
field allowed to take complex values. With this new approach, the dimension, size, 
structure, and periodicity of the lattice play no part at all in the main result. 

Their starting point is that for a real (inverse) temperature, /3, the partition 
function of the finite system is a polynomial in the activity exp(— 2^), where h is 
the external field. Since all the coefficients of the polynomial are positive, none of 
their roots can be on the real positive axis, but are in general complex. 



Lee- Yang Theorem 

The discovery of Lee and Yang is that the zeros of the partition function are all loc- 
ated on the unit circle of the complex activity plane, or equivalently on the imaginary 
/i-axis. This was originally demonstrated for an Ising system with ferromagnetic in- 
teractions (with no need for any limit to first nearest neighbours) although the result 
is valid for more general models |145] . 

The distribution of the zeros on the unit circle will determine whether or not a 
phase transition exists. As the number of spins, N , approaches infinity, if the zeros 
do not condense onto the positive real axis the free energy, F, will remain analytic 
and no phase transition exists at the given /3; otherwise, if the zeros condense onto 
the real positive axis a phase transition will exist at this f3 = Pc- 

We will use the demonstration of the theorem described in |146] including the 
derivation of the magnetisation of the system from the angular distribution of roots 
on the unit circle. We will focus on a general Ising model defined on a graph of 
N sites (vertices of the graphs) with at most one link joining each pair of vertices, 
which are then defined as neighbours. The total number of links is L. 
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The partition function for spins, (Jj, is: 



where is a sum over all links and hi denotes the magnetic site-dependent field 

at the site i (strictly the true external field is proportional to hi/ 13). Let us define 

= 6-2'^' ; T = e-^^. (F.2) 

Then the partition function can be recast in the form 

Zn = ^ exp (^L + P^) , (F.3) 



where P is a polynomial in r and p obtained as 



P = exp 



<Ti=±l 



(F.4) 



For r and pi real and positive, P is always positive and cannot vanish. We assume 
in the following < r < 1, i.e., we are in the ferromagnetic regime. 

It is easy to find that the polynomials P corresponding to the simplest graphs 

are: 

Pi2 = 1 + r(pi + P2) + (P1P2) (F.5) 

P123 = (l+PiT)(l+p3r)+P2(r + pi)(r + p3) (F.6) 



It can then be seen that P is a polynomial of degree one in each pi individually and 
of degree N in all of them. 

We will analyse how the polynomials are generated by building a graph step 
by step. First, one observes that for any disjoint union of subsets of the graph, 
the polynomial P factorizes, P = P1P2, where Pi and P2 are the corresponding 
polynomials of the two separate subsets. With this in hand, we generate a new 
polynomial by identifying a site a from the subset 1 with a site b from subset 2. We 
call P(i2) the corresponding (contracted) polynomial of the union. Since Pi is linear 
in pa and P2 is linear in pj, we can write 

Pi = A^ + PaA_ ; P2 = P+ + pfeP_, (F.7) 

where A^{B^) correspond to the contribution with spin up, cTq = +1 (a?, = +1), 
and A_(P_) correspond to the contribution with spin down, = —1 (o"6 = — 1). 
When the identification of a and b is made, a new activity variable, pab, is attached 
to the new site and we have the following contraction process, see Fig. IF. 11 

P1P2 = A+B+ + PaA_B+ + ptA+B^ + PaPbA.B. (F.8) 
P(i2) = A+B+ + pabA_B_ . (F.9) 
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P2 




Figure F.l: The contraction process obtained by identifying a and h from two disjoint 
subsets of the whole graph. 



Using this process, one can obtain for example the expression of Eq. (1F.6P by 
joining two unconnected two- vertex graphs, such as that of Eq. (IF.Sp . and by identi- 
fying two of the four points. It is a good exercise to obtain the expression for a 
cyclic graph of three vertices by first obtaining that corresponding to a graph of 
four vertices (by joining a graph of two with a graph of three) which is 



^^1234 = l + r{pi[l+p2(l+P3)]+P4[l+P3(l+P2)]} + 

r^[piP4(l + P2 + Ps) + P2(l + Ps) + Ps] + 
T"^(PiP3 + P2P4) +P1P2P3P4, (F-10) 



and then identifying the outer vertices (1 and 4 in the above equation) to obtain for 
the cyclic graph: 



cyclic 



1 + r^(pi + P2 + P3 + P1P2 + P1P3 + P2P3) + P1P2P3 • (F-11) 



The contraction process can be applied to a single connected part, where at first 
the sites a and h are different, then they are identified as a single site ah with activity 
variable pab and 

P = + A^+Pa + A+^Pb + A^^PaPb > Pab = A+++ PabA^^ . (F.12) 

As in the initial set, in the contracted graph no pair of vertices can be joined by 
more than one link. 

Thus we have demonstrated that the contraction process allows the polynomial 
P of any arbitrary graph to be obtained by starting from the elementary result of 
the simplest graph of two vertices joined by a link. In this case, from Eq. flF.51) . the 
zero of the polynomial is obtained for: 

1 + ^P2 

Pi = , , F-13 

T + p2 

which defines a one-to-one mapping between the complex planes pi and p2- For r 
real it leaves the unit circle invariant while for < r < 1 exchanges the interior and 
the exterior of this circle. Therefore it can be stated that if |pi| < 1 and |p2| < 1 , 
or |pi| > 1 and \p2\ > 1, the polynomial can not vanish. 
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The previous property for a graph of two vertices is generahsed in the following: 
For an arbitrary graph, if all pi lie inside, or all pi lie outside the unit circle, P is 
different from zero. 

To demonstrate the foregoing statement, it is sufficient to verify that the property 
survives the contraction process. Let us assume that for a given graph P{pi) 7^ 
when I Pi I < 1 for all i, and apply the contraction process: When the dependence on 
the two points (a and b) to be identified is made explicit, P is given by Eq. (lF.12p . 
while if the identification is already made Pab is a function only of the subset {pi} — 
{Pa, Pb} and a new variable pab- We fix all the p^'s distinct from pa and pb within 
the unit circle. We want to show that \A^^\ > \A |, in which case Pab will be non- 
vanishing for \pab\ < 1. Recalling that P 7^ for \pa\ < 1, \pb\ < 1, the polynomial 

^4++ + p(v4_|__ + v4_+) + p'^A must have its two roots equal to or greater than 

unity, which means that > \A |. Thus Pab is different from zero when all its 

p's are within the unit circle, and this property holds for any graph. 

If we set now all pi = p, corresponding to a uniform external field, from the 
symmetry property under reversal of the field h — )■ —h , p — t- p^^ , and Z{h) = 
Z{—h), we have: 

e+^''P(r, p) = e-^'^P(r, p-^) , (F.14) 

or equivalently: 

PiT,p)=p''Pir,p-'). (F.15) 

As a result if P 7^ for |p| < 1, it follows that also P 7^ for |p| > 1. The Lee- Yang 
theorem has thus been demonstrated - the partition function can (and does) only 
vanish on the unit circle |p| = 1. 

Distribution of Roots on the Unit Circle 

Using the definition of P in Eq. (lF.4p we can obtain the two extreme cases of the 
distribution of the zeros on the unit circle: 

• For a system at infinite temperature (r = 1): 

P(r,p) = P(l,p) = (l + p)^. (F.16) 

• For a system at zero temperature (r = 0): 

P(r,p) = P(0,p) = l + p^. (F.17) 

Then, decreasing the temperature from infinity, one goes from a degenerate zero 
with multiplicity N located at p = —1 to a uniform distribution pk = , In 

addition, if remains finite, no zero will be located on the real positive axis (p = 1, 
h = 0). 

Let us obtain the general result for a lattice of sites, with coordination number 
q, see footnote in Appendix |X1 when taking the infinite volume limit. According to 
Eq. ( ]F.1|) . the free energy per site in a uniform field is: 

F=^AnZ = lq(3 + h + \n2+ lim llnP^(r,p), (F.18) 
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where we have used that L = Nq/2, with L being the total number of hnks in the 
lattice. Pn can be factorized in terms of its roots 



N 



a=l 



n ( 1 - ^ ) ■ 



As was demonstrated in the previous section, the zeros will accumulate for N ^ oo 
on the unit circle, p = e"^, with a (r-dependent) angular density Nfi{ip) satisfying 

/i(<^) = /i(-v9) > ; r d^fiiip) = 1 , (F.20) 

J —n 

where the above property is a consequence of the invariance under field reversal, 
h — )■ —h. Making the change 

^ In P^(r,p) = l|j In (^1--^) ^1 ^ d^Nf,{^)\n(l ^ 



(F.21) 

and using the symmetry property of the angular distribution to join the contributions 
of the conjugate zeros we arrive at 

F = lqf3 + h + \n2 + l- [ d^fi{ip) ln(l + -2p cos if) , (F.22) 

which is valid over the whole range < r < 1. Below the critical temperature, where 
the support of p{ip) is the full circle, this defines in general two different analytic 
functions, one for |p| < 1 and another for |p| > 1. At p = 1 {h = 0), F is continuous. 
Nevertheless, let us consider the behaviour of its derivative with respect to h, i.e. 
the magnetisation, as the temperature varies. By definition 

M = — = l + j d^p{^)— [HI - pe-^n] , (F.23) 

but 

d dp d d 

dh^dhd'p^ %' ^ ^ ^ 

thus 

M = 1 + 2/ dy^M^)-/— = / dM^) ^ (F.25) 

dy^/i(^) . ^ (F.26) 

1 — 2pcosv9 + 

where we made use of the normalisation of p{ip), Eq. (lF.20p . The magnetisation 
vanishes as /i — )■ (p — )■ 1) for r > Tc and becomes discontinuous for r < Tc. The 
above equation is not zero for p = I only if cos<f = 1. Therefore the value of the 
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magnetisation is directly related to the vahic of the density of zeros on the real 
positive axis ((^ = 0), as a simple contour integration summing the residues shows: 

r < Tc ; M± = lim M = ±27r//(0) (F.27) 

We have therefore demonstrated that in the Thermodynamic Limit the spontaneous 
magnetisation is directly related to the existence of zeros on the real positive axis. 
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Appendix G 
IBERCIVIS 



During the year 2007, the BIFI (Institute for Biocomputation and Physics of Com- 
plex Systems) and the National Fusion Laboratory of the CIEMAT (Centro de 
Invest igaciones Energeticas, Medioambientales y Tecnologicas) , collaborating with 
the city council of Zaragoza (Spain), leadered the ZIVIS project. The scope of the 
project was to develop a volunteer supercomputing platform, based on individual 
computers located in both private homes and public buildings, to be used by the 
scientific community in the University of Zaragoza. This network of individual com- 
puters would make it possible to perform calculations as a single installation. The 
project converted Zaragoza into the first city with a large and stable computing 
structure based only on the volunteer effort of its citizens in a non-profit contribu- 
tion to science and research. The scientific goal of the project was the analysis of 
plasma trajectories in the next-generation nuclear fusion reactors. 

The ZIVIS project was based on the BOINC (Berkeley Open Infrastructure for 
Network Computing) software. BOINC was originally developed in the SETI@home 
project for analysing the electromagnetic radiation received from outer space in the 
search for extra-terrestrial intelligence. The software every volunteer downloaded has 
the form of a simple-to-instal program (client) that works as a screensaver during 
idle times of the computer (on average, some 80% of total CPU time). This means 
that the user does not notice any inconvenience when using the computer. A server 
sends tasks via Internet to the clients, which return the results of the computations 
as they are performed and when their Internet connection is enabled. 

ZIVIS was an impressive success that even took its promoters aback. Around 
3000 volunteers offered more than 5000 computers resulting in around 850 000 hours 
of CPU time. It allowed the analysis of more than 4.2 x 10^ ion trajectories. During 
the short life of the project (40 days), it achieved a performance comparable with a 
large "traditional" supercomputing facility (in fact it became one of the five most 
important supercomputing facilities in Spain). 

After the marvellous experience of ZIVIS, and the interest it raised in the rest of 
the country, an extension of the project to a higher level was designed. The result 
was called IBERCIVIS |147] . It was predicted that it would include more than 
50 000 nodes (only in Spain) resulting in the largest computer of this kind all around 
the world. The project was officially launched in June 2008 and it has helped the 
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scientific community since then. The project is becoming a major achievement both 
scientifically (now six applications are producing high-precision numerical results) 
and socially (more than 6000 users are sharing their computers and getting involved 
with scientific research daily). It differs from ZIVIS in some major points: 

• There is not just one scientific application running within IBERCIVIS. Sci- 
entists all around the world are invited to use the infrastructure created to 
run their programs. The "machine" is an open structure where any research 
group could in principle execute their programs. The diagram of the process 
an application must follow to run in IBERCIVIS is outlined in Fig. ICll While 
at first IBERCIVIS started just running three different applications (plasma 
trajectories in fusion reactors, protein folding, and phase transitions in dis- 
ordered systems), it is today running simultaneously six scientific programs 
(the former three plus neuronal amino acid simulations, adsorption in porous 
media, and light behaviour at the nanoscale) [147] . Every volunteer can choose 
which application they prefer to run on their computer. 



Step1 ; selection step 2; porting step 3: production 




Figure G.l: An IBERCIVIS application road map from the beginning of the collab- 
oration up to the dissemination of the results. 



• IBERCIVIS is not a temporary project (unlike ZIVIS or many other BOINC- 
based volunteer computing projects), so it will be possible to submit applica- 
tions indefinitely. This implies that IBERCIVIS is seen by scientists as a stable 
structure like traditional supercomputing centres on which they can run their 
programs via a user-friendly interface with a typical queueing system, launch 
their simulations from their personal computers, and receive the output on a 
special server with high storage capability. 

Advantages and Disadvantages versus Other Struc- 
tures 

On the one hand, the strong points of the project are the following: 

• Apart from the huge scientific interest that nowadays every supercomputing 
facility produces, this one has the additional feature of its extremely low cost. 
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As it mainly uses existing infrastructure (both computers and networks), it 
only requires the effort of the development and support of the specific software 
and servers, apart from an effort in publicity to persuade people to join the 
project (for example, there are periodically competitions between individual 
clients or between teams of clients with prizes for the winners). 

• It provides an excellent way for bringing science close to people; the best 
way to make someone interested in something is to get them involved with it. 
People feel themselves to be part of the solution of a hard research problem 
and learn about the subject. Channels of communication can be established 
between volunteers and scientists through blogs and social networks |147] . In 
this way, the most advanced scientific knowledge is spread to society using 
modern information technologies. 

On the other hand, the main objections to the project, and generally to any kind 
of distributed volunteer computation, are the following: 

• Every node of the supercomputing facility is located, in principle, away from 
other nodes making the communication between nodes very expensive. This 
fact makes direct parallelisation of the computing problem impossible if com- 
munication between nodes is a must. The range of applicability is therefore 
lower than for "one-room" supercomputers. Nevertheless, there exists a large 
class of problem for which communication between nodes is irrelevant. These 
are the so-called embarrassingly parallel problems. Within this class, problems 
can be split into independent simulations whose outputs can be joined later. 

• Given that all the data necessary for the job must be transferred to the volun- 
teer's computer via the Internet, the input/output of the program cannot be 
too massive. Otherwise, it would interfere with the volunteers' network traffic 
and they would naturally become upset. Typically the size of the transfered 
I/O files should not be greater than a few megabytes. 

• Volunteers' computers are not as stable as computers within a traditional 
supercomputing facility - they are more likely to be restarted or even turned 
off. Therefore, in order to increase the probability for a task to be finished, a 
long computation (say of a week) must be divided into short portions (say of 
one hour). This produces an increase both in the network traffic and in the 
probability of corruption of transferred data as they undergo several iterations 
from computer to computer. The implementation of a bulletproof system for 
the validation of the output of each task becomes crucial in order to ensure 
reliability of the final data. 

The Numbers of the Project 

As was said before, the project is being a success both in scientific results and in 
citizen contributions. The number of users and of available nodes are both increasing 
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continuously, see Fig IG.2I The number of registered users is today (May 2009) 
around 11000 and every day around 5500 of them share their computers with the 
scientists. The cluster equivalence of the IBERCIVIS structure is currently of 900 
nodes located in a classical supercomputing facility. 
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Figure G.2: Evolution in time (from June 2008 to June 2009) of the number of 
volunteers, computers, and cores involved in the project. 

Only during the project's first year it has produced around eight millions hours 
of CPU time. The total economic cost has been around 270 000 euros, which is 
not too much taking into account that the first year of life will be surely the most 
expensive. 

At present, six research groups are running simultaneously on the platform, and 
three new applications are in the porting process. This is indicative of the interest 
that the project is producing in the scientific community. 

The publicity campaign of the project has also been really important. Apart from 
appearances in newspapers, magazines, and TV |147] , more than 200 000 entries can 
be found using the Google search engine for the word "ibercivis" . 

Our Experience 

In our case we have been using IBERCIVIS for approximately a year to simulate 
disordered magnetic systems defined in lattices through Monte Carlo methods. In 
particular, we have studied the three-dimensional eight-state Potts model [82j in the 
presence of dilution. The results of these simulations are presented in Sec. 13.4.31 

We started running on IBERCIVIS from its earliest stages (around May 2008) 
so that we have followed all the evolutional process of the project, thus suffering its 
teething problems but also experiencing its gratifying educational advances. 

Our starting point was code written in C that we had been running during the 
previous year on different supercomputing facilities. The code was neither too com- 
plex nor did it make use of any exotic C libraries - facts which made the porting 
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process simpler, ensuring compatibility between the different platforms of the vo- 
lunteer's computers (Windows, Mac, or Linux operating systems with 32 or 64 bit 
processor architectures). 

Our application is the perfect example of an embarrassingly parallelisable one. 
We parallel in four fully independent ways: Firstly we have to simulate different 
system sizes. Secondly, for each size, we have to simulate different values of the 
dilution of the system. Thirdly, for each dilution, we have to simulate different 
realizations of the random spatial hole distribution (each one is called a sample). 
And fourthly, we have to simulate each sample at different values of its internal 
energy. 

Our application does not have strong requirements of RAM within the volunteer 
nodes (around 40 megabytes for the most demanding case) or of disk storage (around 
2 megabytes to store the 1/0 configurations of the largest systems). The main 
problem of our simulations is that, as the system size is increased, the run time 
grows exponentially. This fact forced us to design a continuity system allowing the 
division of long simulations into small (in terms of time) parts. In particular, the 
process that one of our runs for a given dilution and system size follows is: 

1. For each sample (typically there will be around 1000 of them) , a random spatial 
hole configuration is generated depending on the system dilution. The holes 
will remain fixed in time (quenched disorder). 

2. For each energy (typically there will be around 30 of them), to each non- 
empty site of the lattice a spin variable is assigned. The assignation can be 
or a random Potts state, or a fixed one, or even a value depending on the 
system's energy. The result is called a configuration and saved into a file. 

3. Each configuration is sent to a volunteer's machine and updated there using a 
MC method. In addition, some measurements are taken into the system during 
the update process and saved into a measurement file. It has been calculated 
that the optimal time for a run on each volunteer machine (in order to minimise 
both the errors due to unexpected shutdowns and the web traffic) is around 
half an hour. Therefore the number of MC updates of the configuration must 
be set to last around this amount of time. 

4. When the specified mimbcr of MC updates have been made, both the out- 
put configuration and the measurement file are uploaded from the volunteer's 
computer to a server, where they are checked properly in order to avoid cor- 
ruption. 

5. If the total desired simulation time is longer than half an hour, the continu- 
ity system takes over control: It will resend both the configuration and the 
measurement file to another volunteer (step 3), who will continue updating 

the system from the last configuration. The process is repeated until the 
total number of desired MC updates is performed. The number of continuity 
iterations of a given job must be specified a priori in our application. 
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For example, for the smallest simulated system (with 24^ sites) the configurations 
travel only twice between the IBERCIVIS server's and the volunteer's computers 
while for the largest system (with 128'^ sites) the continuity process is repeated 
twenty times. See Fig. IG.3I for a full sketch of the process. 




VOLUNTEERS 



Figure G.3: Different stages of a simulation from the user submission to the results 
download. The term "work unit" is used the define each small portion of the whole 
job that is sent to a single client (each work unit should last around half an hour). 
Both validation and assimilation stages are crucial in order to avoid data corruption. 
The continuity stage is only used if the required time for the job is much longer than 
half an hour. 



The implementation of the continuity systems is by no means naive. Each file 
must be univocally identified in order to avoid misdirections. If only a single file is 
lost or misplaced, the rest of the continuity process will fail for this configuration 
making the analysis of the corresponding sample impossible. In addition, if there 
were a temporary failure in the servers, some files would surely fail to be transferred, 
resulting in a breakdown of the continuity process. This sometimes happened in the 
early stages of the project: due to some server crashes, the continuity process was 
unstable and simulation of the largest systems was impossible. Finally, by building 
paranoid assimilation- validation systems and by developing a univocal nomenclature 
for each file we decreased the failure rate of the continuity system to less than 0.5%. 

Using IBERCIVIS, we have simulated the site-diluted three-dimensional Potts 
model with a precision never reached before. We obtained more than 2.5 x 10^ hours 
of CPU time. A more detailed description of our usage of the infrastructure is given 
in Table EH 
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L 


# dilutions 


# samples 


# iterations 


# energies 


CPU time (xlO^ hours) 


24 


14 


7000 


2 


40 


280 


32 


13 


6500 


2 


30 


195 


48 


12 


24000 


2 


30 


720 


64 


8 


10000 


3 


30 


450 


128 


4 


2000 


20 


30 


900 



Table G.l: Approximate statistics of our simulations in IBERCIVIS for each hnear 
lattice size L. The fourth column gives the number of continuity iterations used to 
obtain enough MC time for each simulated energy. 
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